Arquivo

Posts Tagged ‘normalidade’

Como fazer e interpretar o gráfico quantil-quantil

Gráfico quantil-quatil da distribuição normal com envelopes de confiança simulados.

O gráfico quantil-quantil (q-q) é uma ferramenta muito útil para checar adequação de distribuição de frequência dos dados à uma distribuição de probabilidades. Situações como essa ocorrem principalmente na análise de resíduos de modelos de regressão onde o gráfico q-q é usado para verificar se os resíduos apresentam distribuição normal. O gráfico q-q é melhor que o histograma e o gráfico de distribuição acumulada empírica porque nós temos mais habilidade para verificar se uma reta se ajusta aos pontos do que se uma curva de densidade se ajusta a um histograma ou uma curva de probabilidade acumulada se ajusta à acumulada empírica. Compare às três visualizações com o código a seguir.

#-----------------------------------------------------------------------------
# q-q vs histograma vs ecdf

y <- rnorm(50)
par(mfrow=c(1,3))
qqnorm(y); qqline(y)
plot(density(y))
curve(dnorm(x, mean(y), sd(y)), add=TRUE, col=2)
plot(ecdf(y))
curve(pnorm(x, mean(y), sd(y)), add=TRUE, col=2)

# qual você sente mais segurança para verificar adequação?
# nossos olhos têm mais habilidade para comparar alinhamentos
#-----------------------------------------------------------------------------

Apesar de muito usado, poucos usuários conhecem o procedimento para fazê-lo e interpretá-lo. O procedimento é simples e pode ser estendido para outras distribuições de probabilidade, não apenas para a distribuição normal como muitos podem pensar. Além do mais, alguns padrões do gráfico q-q obdecem à certas características dados dados, como assimetria, curtose e discreticidade. Saber identificar essas características é fundamental para indicar uma transformação aos dados. Abaixo o código para o gráfico q-q para distribuição normal, q-q para qualquer distribuição e o q-q com envelope de confiança obtido por simulação. A execução e estudo do código esclarece o procedimento. O envelope de confiança por simulação torna-se proibitivo para grandes amostras pelo tempo gasto na simulação.

#-----------------------------------------------------------------------------
# função para fazer o gráfico quantil-quantil da normal

qqn <- function(x, ref.line=TRUE){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- qnorm(pteo)           # quantis teóricos sob a normal padrão
  plot(xo~qteo)                 # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- qnorm(c(1,3)/4)       # 1º e 3º quartis teóricos
    points(qrtt, qrto, pch=3)     # quartis, passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
}

x <- rnorm(20)
par(mfrow=c(1,2))
qqn(x)
qqnorm(x); qqline(x)
layout(1)

#-----------------------------------------------------------------------------
# função para fazer o gráfico quantil-quantil de qualquer distribuição

qqq <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1)){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- do.call(distr,        # quantis teóricos sob a distribuição
                  c(list(p=pteo), param))
  plot(xo~qteo)                 # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- do.call(distr,        # 1º e 3º quartis teóricos
                    c(list(p=c(1,3)/4), param))
    points(qrtt, qrto, pch=3)     # quartis, por eles passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
}

x <- rnorm(20)
par(mfrow=c(1,2))
qqq(x)
qqnorm(x); qqline(x)
layout(1)

x <- runif(20)
qqq(x, ref.line=TRUE, distr=qunif, param=list(min=0, max=1))

x <- rgamma(20, shape=4, rate=1/2)
qqq(x, ref.line=TRUE, distr=qgamma, param=list(shape=4, rate=1/2))

#-----------------------------------------------------------------------------
# envelope para o gráfico de quantis (simulated bands)

qqqsb <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1),
                  sb=TRUE, nsim=500, alpha=0.95, ...){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- do.call(distr,        # quantis teóricos sob a distribuição
                  c(list(p=pteo), param))
  plot(xo~qteo, ...)            # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- do.call(distr,        # 1º e 3º quartis teóricos
                    c(list(p=c(1,3)/4), param))
    points(qrtt, qrto, pch=3)     # quartis, passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
  if(sb){
    rdistr <- sub("q", "r",       # função que gera números aleatórios
                  deparse(substitute(distr)))
    aa <- replicate(nsim,         # amostra da distribuição de referência
                    sort(do.call(rdistr, c(list(n=n), param))))
    lim <- apply(aa, 1,           # limites das bandas 100*alpha%
                 quantile, probs=c((1-alpha)/2,(alpha+1)/2))
    matlines(qteo, t(lim),        # coloca as bandas do envelope simulado
             lty=2, col=1)
  }
}

x <- rnorm(20)

#png("f047.png", 400, 300)
qqqsb(x, xlab="Quantis teóricos da distribuição normal padrão",
      ylab="Quantis observados da amostra",
      main="Gráfico quantil-quantil da normal")
#dev.off()

x <- rpois(50, lambda=20)
qqqsb(x, distr=qpois, param=list(lambda=20))

#-----------------------------------------------------------------------------

A normalidade dos resíduos é um dos pressupostos da análise de regressão. São os resíduos e não os dados que devem apresentar normalidade. Se a distribuição dos dados, ou melhor, da sua variável resposta (Y) condicional ao efeito das suas variáveis explicativas for normal, os resíduos terão distribuição normal. Porém, se você aplica um teste de normalidade aos dados (Y) você não está considerando os efeitos das variáveis explicativas, ou seja, você está aplicando um teste na distribuição marginal de Y que não tem porque atender a normalidade. Todo teste de normalidade supõe que os dados têm uma média e uma variância e só os resíduos atendem essa premissa porque os dados (Y) têm média dependente do efeito das variáveis explicativas.

Distribuição condicional e marginal da variável resposta.

#-----------------------------------------------------------------------------
# distribuição condicional vs distribuição marginal

layout(1)
x <- rep(1:4, e=50)
y <- rnorm(x, mean=0+0.75*x, sd=0.1)
da <- data.frame(x, y)
plot(y~x, da)

db <- split(da, f=da$x)

#png("f046.png", 400, 300); par(mar=c(1,1,1,1))
plot(y~x, da, xlim=c(1,6), xaxt="n", yaxt="n", xlab="", ylab="")
abline(a=0, b=0.75)
lapply(db,
       function(d){
         dnst <- density(d$y)
         lines(d$x[1]+dnst$y*0.1, dnst$x)
         abline(v=d$x[1], lty=2)
       })
points(rep(5, length(y)), da$y, col=2)
abline(v=5, lty=2, col=2)
dnst <- density(da$y)
lines(5+dnst$y*2, dnst$x, col=2)
text(3, 1, "distribuição\ncondicional")
text(5.5, 1, "distribuição\nmarginal", col=2)
#dev.off()

par(mfrow=c(1,2))
qqnorm(da$y)
qqnorm(residuals(lm(y~x, da)))
layout(1)

#-----------------------------------------------------------------------------

Muitos usuários preferem aplicar um teste de normalidade do que olhar para o gráfico q-q. Isso têm duas razões: (1) costume, o usuário sempre usou aplicativos para análise de dados que não dispõem de recursos gráficos, eles conduzem toda análise sem sequer ver os dados em um gráfico, (2) consideram subjetiva a análise gráfica. Do meu ponto de vista, a subjetividade está presente ao aplicar um teste também pois o usuário é quem escolhe o teste e o nível de significância. Outro ponto é que os testes de normalidade assumem independência entre as observações e os resíduos não são independentes e foram obtidos após a estimação de p parâmetros, o que não é considerado por nenhum teste. Ou seja, qualquer teste aplicado fornece um resultado aproximado. Mas o que de fato eu defendo é que a análise gráfica é indiscutivelmente mais informativa. Veja, se o teste rejeita a normalidade é porque os dados não apresentam distribuição normal por algum motivo. Quando você visualiza o q-q é possível explicar a fuga de normalidade que pode ser sistemática: (1) devido à desvio de assimetria, (2) de curtose, (3) à mistura de distribuições, (4) à presença de um outlier e (5) ao dados serem discretos. Essas são as principais causas de afastamento. Cada uma delas sugere uma alternativa para corrigir a fuga: tranformação, modelagem da variância, deleção de outlier, etc. Nesse sentido, como identificar esses padrões de fuga? É o que os gráficos animados vão mostrar. Até a próxima ridícula.

Fuga de normalidade por assimetria representada pelo gráfico quantil-quantil. A assimetria no q-q aparece com pontos dispostos em forma de arco.

#-----------------------------------------------------------------------------
# assimetria

dir.create("frames")
setwd("frames")

png(file="assimet%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in 10*sin(seq(0.01, pi-0.01, l=100))){
  curve(dbeta(x, i, 10-i), 0, 1, xaxt="n", yaxt="n", xlab="", ylab="")
  y <- rbeta(100, i, 10-i)
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 assimet*.png assimet.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

#-----------------------------------------------------------------------------

Fuga de normalidade por mistura de distribuições (alteração de curtose) representada pelo gráfico quantil-quantil. Mistura de distribuições geram fugas nos extremos.

#-----------------------------------------------------------------------------
# mistura de distribuições

png(file="mistura%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in sin(seq(0, pi, l=100))){
  curve(i*dnorm(x,0,1)+(1-i)*dnorm(x,0,6), -20, 20,
        xaxt="n", yaxt="n", xlab="", ylab="")
  y <- c(rnorm(ceiling(500*i),0,1), rnorm(floor(500*(1-i)),0,6))
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 mistura*.png mistura.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

#-----------------------------------------------------------------------------

Fuga de normalidade devido dados serem discretos (padrão escada) representada pelo gráfico quantil-quantil.

#-----------------------------------------------------------------------------
# discreticidade

png(file="discret%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in c(1:100, 99:1)){
  x <- 0:i
  px <- dbinom(x, i, 0.5)
  plot(px~x, type="h", xaxt="n", yaxt="n", xlab="", ylab="")
  y <- rbinom(100, i, 0.5)
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 discret*.png discret.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

#-----------------------------------------------------------------------------
Anúncios

Análise de resíduos em regressão não linear

Conjunto de gráficos para análise dos resíduos de um modelo de regressão não linear.

A análise dos resíduos de um modelo é feita para verificar a plausividade das pressuposições envolvidas. Os modelos linerares de regressão clássico, ou seja, aqueles em que as observações são realizações independentes (independência) e apresentam a mesma dispersão (homogeneidade de variância), podem ser ajustados no R com a função lm() e aov(). Para objetos dessas classes existe um método plot.lm() que apresenta os gráficos de análise de resíduo. Porém, modelos de regressão não linear podem ser ajustados com a função nls() que não possui um método para análise de resíduos.

Para aplicação de inferência (testes de hipótese, intervalos de confiança, etc), o modelo não linear é aproximado linearmente. Com isso eu quero dizer que podemos usar as mesmas técnicas de análise de resíduos para modelos de regressão não linear. A matriz gradiente do modelo (de derivadas da função em relação ao vetor de parâmetros) pode ser usada dentro da função lm() e com isso podemos obter facilmente os gráficos de análise de resíduos.

Esse procedimento ainda permite obter a análise de variância para o modelo de regressão não linear, fazer a partição ortogonal da soma de quadrados total em devido ao modelo de regressão e devido aos desvios de regressão. No entanto, esse não deveria ser o quadro apresentado se compararmos o que obtemos com modelo de regressão linear. No modelo linear particionamos a soma de quadrados total corrigida para o intercepto. O nosso modelo de regressão não linear se torna um modelo de intercepto se V e D forem zero (ver modelo do CMR). Portanto, o modelo nulo seria aquele em que estimamos apenas A e que podemos comparar com o modelo completo por meio da função/método anova.nls(). Esse é o teste de hipótese para o modelo de regressão não linear. Vale lembrar que alguns modelos não lineares não podem ser reduzidos a modelos de intercepto via restrição paramétrica. Para esses modelos deve-se ter cuidados ao usar o R², inclusive.

Abaixo apresento o procedimento para ajuste de um modelo não linear, extração da matriz gradiente e obtenção dos gráficos de resíduos. No final eu fiz uma função (R2()) para facilitar a obtenção do quadro de análise de variância e R². Não é porque você tem uma função para isso que vai sair aplicando ela em qualquer modelo não linear, ok? Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de postássio liberado em função do tempo

klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
                     91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
                     89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
                     104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
                     128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
                     132.67, 154.48, 129.11, 151.83, 147.66, 127.30),
                   t=rep(c(15, 30, 45, 60, 75, 90,
                     120, 150, 180, 210, 240, 270), each=3))

#-----------------------------------------------------------------------------
# ajustando o modelo de regressão não linear aos dados

n0 <- nls(k~A*t/(V+t)+D*t, data=klib, start=list(A=90, V=15, D=0.21))
summary(n0)

#-----------------------------------------------------------------------------
# extraindo a matriz gradiente avaliada nas estimativas dos parâmetros

F <- attr(n0$m$fitted(), "gradient")
F

#-----------------------------------------------------------------------------
# passando a matriz gradiente para a lm(), importante remover intercepto

m0 <- lm(k~-1+F, data=klib)

#-----------------------------------------------------------------------------
# gráfico de análise dos resíduos

#png("f007.png", w=500, h=500);
par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1))
plot(m0)
mtext("Análise de resíduos para modelo de regressão não linear",
      outer=TRUE, line=-2, cex=1.4) 
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# veja que o ajuste é o mesmo pelas medidas abaixo

cbind(fitted(m0), fitted(n0))       # valores ajustados
cbind(residuals(m0), residuals(n0)) # valores preditos
c(summary(m0)$sig, summary(n0)$sig) # estimativa do desvio padrão residual
vcov(m0); vcov(n0)                  # matriz de covariância das estimativas

#-----------------------------------------------------------------------------
# quadro de anova com SQ de regressão e SQ de resíduo

anova(m0) # partição da soma de quadrados total
anova(n0, lm(k~1, klib)) # SQ do modelo não linear corrigido para intercepto

#-----------------------------------------------------------------------------
# função que retorna a anova e R2 para modelos de regressão não linear

R2 <- function(nls.obj){
  da <- eval(nls.obj$data)
  resp.name <- all.vars(summary(nls.obj)$formula)[1]
  form <- paste(resp.name, "~1", sep="")
  m0 <- lm(form, da)
  an <- anova(nls.obj, m0)
  sqn <- deviance(nls.obj)
  sqe <- deviance(m0)
  r2 <- 1-(sqn/sqe)
  aov <- data.frame(fv=c("regression","residuals"),
                    gl=c(-an$Df[2],an$Res.Df[1]),
                    sq=c(-an$Sum[2],an$Res.Sum[1]))
  aov$qm <- aov$sq/aov$gl
  aov$F <- c(aov$qm[1]/aov$qm[2], NA)
  aov$"Pr(>F)" <- c(1-pf(aov$F[1], df1=aov$gl[1], df2=aov$gl[2]), NA)
  names(aov) <- c(" ","Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
  return(list(anova=aov, R2=r2))
}

R2(n0)

#-----------------------------------------------------------------------------