Arquivo

Archive for the ‘inferência’ Category

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

Interpretação da matriz de covariância das estimativas

Quero falar de um erro comum de interpretação de resultados em análise de regressão através de um exemplo. Considere que você tem pessoas com diferentes pesos e diferentes alturas. Facilmente você aceita que uma pessoa mais alta tem maior peso, certo? Ou seja, existe uma correlação positiva entre peso e altura. Pois bem, vamos simular observações desse experimento e em seguida ajustar uma regressão linear simples para relação peso-altura.

require(MASS)

set.seed(12345)
da <- as.data.frame(mvrnorm(80, c(altura=30, peso=60),
                            matrix(c(2,1.9,1.9,3), 2,2)))
plot(peso~altura, da)

m0 <- lm(peso~altura, data=da) # ajusta a reta
summary(m0)              # estimativa dos parâmetros
abline(m0)               # adiciona uma reta ao gráfico
vcov(m0)                 # covariância das estimativas

Aqui está o ponto que eu quero comentar: a interpretação da matriz de covariância. Perceba que a covariância foi negativa. Já vi gente interpretando isso da seguinte forma: espera-se que pessoas mais altas tenham menor peso. Duplamente errado. Primeiro nos sabemos por experiência que a correlação de peso e altura é positiva. Segundo, a matriz de covariância se refere as estimativas dos parâmetros e não as variáveis envolvidas. Nunca a interprete dessa forma.

Então como interpretar? Bem, a matriz de covariância das estimativas, é um reflexo da função objetivo ao redor da solução. A função objetivo nesse caso é minimizar a soma de quadrados (SQ). Então se eu aumento o valor do parâmetro b0, o parâmetro b1 tem que diminuir para compensar, ou seja, para diminuir a SQ. Ocorre um efeito compensatório aqui. Além do mais, esse efeito pode ser eliminado facilmente aplicando uma reparametrização, como por exemplo, centrar a covariável na média.

da$alturac <- with(da, altura-mean(altura)) # centrando a covariável
par(mfrow=c(1,2))
plot(peso~altura, da)
abline(m0)
plot(peso~alturac, da)

m1 <- lm(peso~alturac, data=da) # ajusta a reta
summary(m1)              # estimativa dos parâmetros
abline(m1)               # adiciona uma reta ao gráfico
vcov(m1)                 # covariância das estimativas

require(ellipse)
par(mfrow=c(1,2)) # regiões de confiança conjunta
plot(ellipse(vcov(m0)), type="l")
plot(ellipse(vcov(m1)), type="l")

cov2cor(vcov(m0))
cov2cor(vcov(m1))

Perceba que os resultados são os mesmos em termos de estatísticas de teste e medidas de ajuste. Isso é esperado por eu só fiz uma traslação da altura. Mas o importante é que agora a matriz de covariância tem covariâncias praticamente nulas, que é resultado da translação. O intercepto é então o valor esperado para alturac igual a zero (centro dos dados). Veja porque essa covariância é zero: se eu aumentar b0, não adiante eu alterar b1 que não vai diminuir a SQ, e vice-versa, porque agora eles são ortogonais. Ortogonalidade entre parâmetros é uma característica desejável pois permite você inferir sobre um deles sem considerar o outro. Além disso, tem vantagens do ponto de vista de estimação por métodos numéricos de otimização. Em outras palavras, pegando conceitos de probabilidade, se a distribuição amostral de duas variáveis aleatórias é normal bivariada com covariância nula, a distribuição condicional de A|B é igual a distribuição marginal de A pela independência. Reforçando, eu não preciso conhecer valores de B para descrever A. Considerando tudo que foi argumentado, é sempre preferível que você adote o modelo que apresente menor covariância entre os parâmetros. Até a próxima ridícula.

Categorias:inferência Tags:, ,

Gamma não central (ncg v 0.1.0)

A distribuição gama incompleta não-central pode ser vista como a forma geral da distribuição \chi^2 não-central e é expressa como uma mistura de uma Poisson com a função gamma (\Gamma) incompleta.

Alguns autores propuseram aproximações para a obtenção da função gamma incompleta não-central, sem, no entanto, que os procedimentos numéricos tivessem sido implementados. Outros trabalhos apresentam métodos para a obtenção da função de distribuição gamma (\Gamma) incompleta não-central e de sua inversa, mas os códigos de suas soluções não estão disponíveis.

Para o caso particular da distribuição \chi^2 não-central, existem soluções numéricas de alta qualidade e implementações em programas como o R (família de funções *chisq(), usando-se o argumento ncp). Entretanto, um dos problemas ainda remanescentes é o da obtenção do parâmetro de não-centralidade, dados os valores da função de distribuição, de x e de \alpha, tanto para a distribuição gama (\Gamma), quanto para a \chi^2 não-central.

Neste post apresentamos o pacote (ncg) que implementa uma adaptação ao método originalmente proposto para a função \beta incompleta não-central e um outro que combina o método desses autores com o método da inversão da função de distribuição em relação ao parâmetro de não-centralidade, utilizando para isso o método de Newton-Raphson.

Maiores detalhes acerca da teoria e eficácia do implementado no pacote estão disponíveis na dissertação que serviu de cerne ao desenvolvimento do Pacote, de autoria de Izabela Regina Cardoso de Oliveira, sob orientação do Prof Dr. Daniel Furtado Ferreira, do programa de Pós-Graduação em Estatística e Experimentação Agrícola da Universidade Federal de Lavras (texto completo aqui). Detalhes do pacote, estão também disponíveis no CRAN, (acesse-o por aqui)

Para exemplificar, suponha uma v.a. Y com distribuição \chi^2 com \nu = 5 e parâmetro de não-centralidade \delta = 1. Para calcular a probabilidade acumulada em y = 4 usa-se a função pchisq(), como segue:

## install.packages('ncg')
require(ncg)

nu <- 5
y <- 4
delta <- 1

pchisq(y, nu, 1)

Como a \chi^2 não-central é um caso particular da distribuição \Gamma não-central, por meio de uma simples transformação podemos calcular a probabilidade acumulada da \chi^2, usando a função pgammanc() do pacote ncg, fazemos portanto:

alpha <- nu / 2
x <- y / 2

pgammanc(x, alpha, delta) 

A função deltagammanc() computa o parâmetro de não-centralidade da função \Gamma não-central dados os valores de x, \alpha e p. Como a \chi^2 não-central é um caso particular da \Gamma não-central, essa função também pode ser usada para obter o parâmetro de não-centralidade de uma \chi^2 não-central.

deltagammanc(x, alpha, 0.3471897)

Espero que tanto o pacote, quanto as breves explicações aqui ditas sejam úteis. Agradeço aos colegas que estiveram envolvidos, Izabela e Prof Daniel, obrigado!

Fuzzy c-means, um exemplo!

Pois bem, por sugestão do Professor Julio da Motta Singer apresentarei um script para realizar uma análise de agrupamento no escopo da lógica fuzzy. Seja um experimento fatorial (a mesma natureza dos dados aqui disponibilizados) a seguinte implementação realiza a clusterização de uma dado fator marginalmente ao outro fator.

A heurística da lógica fuzzy tem um incremento interpretativo de, ao contrário do agrupamento k-means, que é dicotômico, associar a cada cluster uma “probabilidade” de associação. Todavia, no método arbitramos uma constante (m) de fuzzyficação. Assim o seguinte script realiza o agrupamento fuzzy (por intermédio da função cmeans() da biblioteca e1071) para vários valores m assim é possível observar o comportamento (decaimento) do erro associado a clusterização e também a disposição dos fatores em relação ao cluster.

Infelizmente (ou para alegria geral) não vou me atrever a discorrer da teoria envolvida nessa análise e um estudo específico deixo a cargo dos maiores interessados.

Todo esse mise en cene, me soa um pouco genérico demais, eu sei (–Me desculpem!). Mas acredito que o script possa ser útil de alguma forma.

Aproveito o ensejo para agradecer as sugestões do Prof Julio da Motta Singer ao trabalho e o incentivo a esse post. Também deixo um alô ao colegas que participaram, de uma forma ou de outra!

Vamos ao trabalho, primeiro instalando (para quem ainda não o fez) e carregando a biblioteca e o conjunto de dados!

## install.packages('e1071')
library(e1071) # carrega pacote com implementação 'c-means'

dados <- transform(read.table("http://dl.dropbox.com/u/38195533/fac_exemplo.txt", # arquivo
                              header = TRUE,
                              sep = '\t',
                              na.string = '.',
                              dec = ','),
                   fatorA = factor(fatorA),
                   fatorB = factor(fatorB))

str(dados)

ajuste <- lm(resp ~ fatorA + fatorB, data = dados)
anova(ajuste) # ANOVA

matriz <- tapply(dados$resp, list(dados$fatorB, dados$fatorA), mean)

A primeira etapa foi realizada. Não posso passar ou próximo trecho sem deixa-los a par do quem está por vir! O agrupamento que “tentamos” é com relação ao um padrão de comportamento do fator sob análise, em relação ao outro fator, assim constituímos ideótipos e assim procedemos. O modo como constituimos esses ideótipos é aplicado unicamente a natureza desse conjunto de dados, mas outros ideótipos, ao seu critério podem ser usados, fique a vontade!

Antes de passarmos agora a análise em si manipulamos os dados e constituímos os centróides:

FinWil <- apply(matriz, 2, mean) - mean(matriz) # Finlay & Wilkinson

altos <- which(FinWil >= 0) # ambientes acima da media
baixo <- which(FinWil < 0) # ambientes abaixo da media

medias <- cbind(apply(matriz[, altos], 1, mean),
                apply(matriz[, baixo], 1, mean))

##--------------------------------------------------------------------
## IDEOTIPOS
## I   - maximo em todos 
## II  - maximo nos altos minimo nos baixos
## III - minimo nos altos maximo nos baixos
## IV  - minimo em todos 
## ADICIONAIS:
## V   - media em todos
## VI  - maximo nos altos e media nos baixos
## VII  - media nos altos maximo nos baixos

I <- c(max(matriz[, altos]), max(matriz[, baixo]))
II <- c(max(matriz[, altos]), min(matriz[, baixo]))
III <- c(min(matriz[, altos]), max(matriz[, baixo]))
IV <- c(min(matriz[, altos]), min(matriz[, baixo]))
V <- c(mean(matriz[, altos]), mean(matriz[, baixo]))
VI <- c(max(matriz[, altos]), mean(matriz[, baixo]))
VII <- c(mean(matriz[, altos]), max(matriz[, baixo]))

centroides <- rbind(I, II, III, IV, V, VI, VII) # centróides

Resta agora realizar a análise em si, veja que o laço for() repete a análise para vários valores da constante m e de cada uma, armazena as “probabilidades” de associação a cada centróide e o erro do ajustamento!


ms <- seq(1.2, 10, by = .1)
lista <- array(NA, dim = c(nlevels(dados$fatorB), 
                           nrow(centroides), 
                           length(ms)))
erros <- matrix(NA, nrow = length(ms), ncol = 1)
dimnames(lista) <- list(levels(dados$TRAT), 
                        c('I', 'II', 'III', 'IV', 'V', 'VI', 'VII'), 
                        ms)

for(i in 1:length(ms))
  {
    j <- ms[i]
    ## função de cluster 'fuzzy c-means'
    grupamentos <- cmeans(medias, # matriz dos dados
                          centroides, # centróides dos clusters
                          verbose = FALSE, # não imprime histórico
                          dist = 'euclidean', # distância empregada
                          m = j, # constante de desfuzzyficação
                          method = 'cmeans') # método
    lista[,,i] <- grupamentos$membership
    erros[i,] <- grupamentos$withinerror
  }

Pronto… feito, quer dizer, quase! Realizamos portanto a análise para os vários valores de m resta agora “escolher” um deles e partir para o resto! Até a próxima…

Ilustração do método Newton-Raphson em regressão não linear

Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação de parâmetros em um modelo de regressão não linear.

Em regressão não linear, as estimativas dos parâmetros são obtidas por emprego de algum método numérico. Um dos métodos empregados é o de Newton-Raphson. Tal método é usado para encontrar raízes de uma função. Portanto, não é um método de maximação/mínimização de funções, mas usa o fato equivalente de que a derivada no ponto crítico de uma função é zero.

Documentação sobre o método Newton-Raphson existe em abundância. O que motivou fazer esse post foi não ter encontrado nenhum que exemplificasse aplicação para estimação de parâmetros em um modelo de regressão não linear. O estimador de mínimos quadrados para um modelo de regressão é aquele que torna mínima a função

S(\theta) = \sum_{i=1}^n (y_i-f(x_i, \theta))^2.

O método Newton-Raphson vai procurar as raízes da derivada dessa função com relação à cada parâmetro. Então, se \theta é um vetor de p parâmetros, a derivada de S(\theta) será uma função p-dimensional,

S_1(\theta) = \displaystyle\frac{\partial S(\theta)}{\partial\theta}=\begin{bmatrix}  \displaystyle\frac{\partial S(\theta)}{\partial\theta_1}\\  \vdots\\  \displaystyle\frac{\partial S(\theta)}{\partial\theta_p}  \end{bmatrix}.

O método se baseia na aproximação de primeira ordem da função S_1(\theta) de forma que um valor da função para um valor \theta^{i+1} na redondeza de um valor \theta^i é

S_1(\theta^{i+1}) \approx S_1(\theta^i)+ (\theta^{i+1}-\theta^i)\displaystyle\frac{\partial S_1(\theta)}{\partial\theta}|_{\theta=\theta^i}.

O valor de S_1(\theta^{i+1}) é zero uma vez que \theta^{i+1} é a raíz da função considerando a expansão de primeira ordem, assim, reorganizando a expressão, temos que a atualização dos valores é data por

\theta^{i+1} = \theta^i- \displaystyle\frac{S_1(\theta^i)}{S_2(\theta^i)},

sendo que

S_2(\theta) = \frac{\partial S_1(\theta)}{\partial\theta}=\begin{bmatrix}  \frac{\partial S_1(\theta)}{\partial\theta_1\partial\theta_1} & \ddots\\  \ddots& \frac{\partial S_1(\theta)}{\partial\theta_p\partial\theta_p}  \end{bmatrix}

é a matriz de derivadas que tem dimensão p\times p.

Esse procedimento numérico foi implementado no script abaixo para estimação dos parâmetros do modelo Michaelis-Menten e o modelo Exponencial reparametrizado que desenvolvi na minha Dissertação de Mestrado. No início eu monto as funções separadas e no final eu faço uma função que recebe os ingredientes mínimos e faz a estimação dos parâmetros. No final do script tem um procedimento gráfico para seleção dos chutes iniciais e acompanhar a trajetória durante as iterações. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# método Newton-Raphson
# objetivo é mínimizar a soma de quadrados dos resíduos mas iremos
# encontrar as raízes da derivada dessa função pelo método de Newton-Raphson

#-----------------------------------------------------------------------------
# parâmetros, modelo e dados simulados deles

theta <- c(A=10, B=2) # parâmetros
x <- 1:12
f <- function(theta, x){
  ## função que retorna os valore preditos para theta e x fornecidos
  ## theta: vetor de parâmetros A e B do modelo Michaelis-Menten
  ## x: vetor de valores da variável independente
  haty <- theta[1]*x/(theta[2]+x)
  return(haty)
}
set.seed(222); e <- rnorm(x,0,0.1)
y <- f(theta=theta, x=x)+e
plot(y~x)
curve(f(theta, x), add=TRUE, col=2)

#-----------------------------------------------------------------------------
# função objetivo, derivada primeira e segunda com relação à theta

S0 <- function(theta, f, y, x){
  ## função que retorna a soma de quadrados dos resíduos para
  ## theta: vetor de parâmetros do modelo f
  ## f: função do modelo de regressão não linear que depende de theta e x
  ## y e x: vetores dos valores das variáveis dependente e independente
  sqr <- sum((y-f(theta, x))^2)
  return(sqr)
}

S0(theta=theta, f=f, y=y, x=x)

# derivada primeira, dimensão 1 x p
D(expression((y-A*x/(B+x))^2), "A") # derivada de S0 com relação a A
D(expression((y-A*x/(B+x))^2), "B") # derivada de S0 com relação a B

# derivada segunda, dimensão p x p
D(D(expression((y-A*x/(B+x))^2), "A"), "A") # derivada de S0 com relação a A e A
D(D(expression((y-A*x/(B+x))^2), "B"), "B") # derivada de S0 com relação a B e B
D(D(expression((y-A*x/(B+x))^2), "A"), "B") # derivada de S0 com relação a A e B
D(D(expression((y-A*x/(B+x))^2), "B"), "A") # derivada de S0 com relação a B e A

#-----------------------------------------------------------------------------
# para ter essas funções vamos usar a deriv3()

ff <- deriv3(~(y-A*x/(B+x))^2,
             c("A","B"),
             function(x, y, A, B){ NULL })

ff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# mas as matrizes são 1 x p e p x p, então falta fazer a soma em n

fff <- function(x, y, A, B){
  ## função que retorna gradiente e hessiano avaliado em theta
  ## y e x: vetores dos valores das variáveis dependente e independente
  ## A e B: valores dos parâmetros do modelo Michaelis-Menten
  m <- ff(x, y, A, B)
  g <- attr(m, "gradient"); G <- apply(g, 2, sum)
  h <- attr(m, "hessian"); h <- matrix(h, ncol=2*2)
  H <- apply(h, 2, sum); H <- matrix(H, 2, 2)
  return(list(G=G, H=H))
}

fff(x=x, y=y, theta[1], theta[2])

#-----------------------------------------------------------------------------
# agora é só fazer o Newton-Raphson

th0 <- c(9,2); tol <- 1e-10; niter <- 100;
dis <- 1; i <- 0
while(dis>tol & i<=niter){
  GH <- fff(x=x, y=y, th0[1], th0[2])
  G <- GH[[1]]
  iH <- solve(GH[[2]])
  th1 <- th0-iH%*%G
  dis <- sum((th0-th1)^2)
  i <- i+1
  print(c(i, dis, th1))
  th0 <- th1
}

n0 <- nls(y~A*x/(B+x), start=c(A=9, B=2))
coef(n0)

iH/vcov(n0) # hessiano é proporcional à matriz de covariância dos parâmetros

#-----------------------------------------------------------------------------
# vamos incluir todos os elementos em uma única função!

newton <- function(formula, data, start, tol=1e-10, niter=100, trace=TRUE){
  ## formula: expressão da resposta como função das covariáveis e parâmetros
  ## data: data.frame que contém dos dados
  ## start: vetor ou lista nomeados de valores iniciais para os parâmetros
  ## tol: norma mínima entre valores consecutivos dos parâmetros
  ## niter: número máximo de iterações do algortimo
  ## trace: se TRUE imprime o histórico de iterações
  form <- paste("~(",paste(formula[2:3], collapse="-"),")^2", sep="")
  vars <- all.vars(formula)
  nvars <- length(vars)
  pars <- 1:(nvars+1)
  l <- as.list(pars)
  names(l)[1:nvars] <- vars
  f <- as.function(l)
  body(f) <- NULL
  ff <- deriv3(as.formula(form), names(start), f)
  th0 <- unlist(start)
  i <- 0; dis <- 100; iter <<- matrix(th0, ncol=length(start))
  while(dis>tol & i<=niter){
    m <- do.call(ff, c(as.list(data), as.list(th0)))
    g <- attr(m, "gradient"); G <- apply(g, 2, sum)
    h <- attr(m, "hessian"); h <- matrix(h, ncol=length(start)^2)
    H <- apply(h, 2, sum); H <- matrix(H, nrow=length(start))
    iH <- solve(H)
    th1 <- th0-iH%*%G
    dis <- sum((th0-th1)^2)
    i <- i+1
    iter <<- rbind(iter, t(th1))
    if(trace==TRUE) print(c(i, dis, th1))
    th0 <- th1
  }
  th0 <- as.vector(th0)
  names(th0) <- names(start)
  ## retorna o último valor iterado e a uma matriz com o histórico
  return(list(theta=th0, iter=iter))
}

#-----------------------------------------------------------------------------
# vamos estudar vetor de chutes presentes nessa lista

start.l <- list(c(A=10,B=2.5), c(A=10,B=0.1), c(A=14,B=2.5),
                c(A=15.13, B=1.7), c(A=18.4,B=0.13))

data <- data.frame(x=x, y=y)
formula <- y~A*x/(B+x)
start <- start.l[[4]]
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-3,15,l=100)
B.grid <- seq(-3,15,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
#png("f023.png", w=500, h=500)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#dev.off()

#-----------------------------------------------------------------------------
# com o mouse selecione regiões da superfície para serem valores iniciais

loc <- locator()
loc <- as.data.frame(loc)
names(loc) <- c("A","B")
start.l <- split(loc, f=1:nrow(loc))

#-----------------------------------------------------------------------------
# ver como fica com o modelo Exponencial

n0 <- nls(y~A*(1-exp(-log(2)*x/B)), start=c(A=9, B=2))
coef(n0)

data <- data.frame(x=x, y=y)
formula <- y~A*(1-exp(-log(2)*x/B))
start <- list(A=5, B=0.5)
nr <- newton(formula, data, start)
nr$theta

sqe <- function(A, B, y, x){ hy <- A*(1-exp(-log(2)*x/B)); sum((y-hy)^2) }
SQE <- Vectorize(sqe, c("A", "B"))

A.grid <- seq(-10,20,l=100)
B.grid <- seq(-10,20,l=100)
sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
        xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
for(i in seq(nrow(iter)-1)){
  arrows(iter[i,1], iter[i,2],
         iter[i+1,1], iter[i+1,2],
         col=2, length=0.1)
  Sys.sleep(0.1)
}
#-----------------------------------------------------------------------------

Verossimilhança para um modelo de efeitos aleatórios

Contornos perfilhados de deviance para regiões de 0,9; 0,95; e 0,99 de confiança. O ponto + no gráfico representa os valores paramétricos usados para simulação dos dados. Imperfeições nos contornos apontam problemas numéricos associados ao passo de integração.

Modelos com termos de efeito aleatórios estão sendo cada vez mais usados na modelagem de dados. Nestes modelos, parte da variabilidade observada na resposta pode ser explicada por um conjunto de fatores de efeito fixo, outra parte pode estar associada a realização de variáveis latentes (ou não observáveis)  a qual é atribuída o efeito aleatório.

A decisão sobre o que é fixo ou aleatório nem sempre é bem clara e depende dos objetivos da pesquisa/inferência. Um caso bem típico é a questão dos blocos nos experimentos agronômicos. Uma maneira simples de pensar é questionar: os blocos presentes no meu experimento foram sorteados de um universo de blocos possíveis? Se sim, considere eles como de efeito aleatório. Há casos em que eles não são de efeito aleatório como, por exemplo, quando o bloco identifica o nível topográfico das parcelas, a classe de peso dos animais. Há casos que é genuinamente de efeito aleatório como, por exemplo, o lote da sementes em teste. Mas como eu ressaltei, a decisão sobre fixo/aleatório depende dos objetivos e difere de pessoa para pessoa.

Quando você declara um fator com efeito aleatório, você assume que os efeitos que ele exerce sobre a resposta são realizações de uma distribuição de probabilidades. Imagine um experimento para verificar a germinação de sementes onde temos 10 lotes de sementes e 4 repetições por lote. Os lotes estudados foram sorteados de uma grande lista de lotes. Ao realizar o experimento, verifica-se variabilidade entre repetições de um mesmo lote e entre lotes. Se lote fosse de efeito fixo, teríamos que estimar 10 parâmetros, um para cada lote. Se lote é de efeito aleatório temos apenas que estimar os parâmetros da distribuição de probabilidades associada a realização dos efeitos. Normalmente é usada a distribuição normal, e nesse caso é estimado a variância dos efeitos aleatórios.

Os efeitos aleatórios não são parâmetros, são realizações de uma distribuição de probabilidades. A verossimilhança para um modelo com efeitos aleatórios é baseado nessa afirmativa. A distribuição marginal da resposta é a integral no domínio dos efeitos aleatórios do produto entre a distribuição condicional da resposta ao valor realizado dos efeitos e a densidade dos efeitos aleatórios, ou seja

f(y) = \int f(y|b) \times f(b)\,\, \text{d}b

em que y é a resposta observada e b o efeito aleatório. Para o nosso experimento com germinação de sementes, temos uma resposta binomial (sementes germinadas em ensaios de 50 sementes), 10 níveis de lote e 4 repetições por lote. A verossimilhança seria

f(\sigma^2_l;y) = \prod_i \int \prod_j \text{Binomial}(y_{ij}|b_i) \times \text{Normal}(b_i, \sigma^2_l)\,\, \text{d}b_i.

No script abaixo eu fiz a simulação de dados para esse modelo, a função de verossimilhança e a otimização usando a optim() bem como a obtenção das regiões de confiança baseados na deviance perfilhada. As imperfeições observadas nos contornos de deviance são causadas por problemas numéricos na etapa de integração. Para checagem fiz o ajuste pela função lme4::glmer(), porém, esta função usa aproximação de Laplace (por isso estimativas levemente diferentes) e possivelmente usa uma função de verossimilhança sem as constantes padronizadoras (por isso máximo da verossimilhança tão diferente).

Para o caso de dados com distribuição Normal, o passo de integração é evitado porque a distribuição marginal possui forma fechada. Além do mais, diversos métodos numéricos podem ser empregados no passo de integração, como aproximação de Laplace, quadratura Gaussiana, integração Monte Carlo.  Outra coisa que é importante são as parametrizações usadas. Em geral, recomenda-se que os parâmetros tenham domínio nos reais. Para o caso da variância é usual otimizar o log do desvio-padrão.

Nos próximos posts vou abordar o perfil de verossimilhança para um parâmetro e a predição dos efeitos aleatórios. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados gerados de um experimento de germinação de sementes

sigma <- 1 # desvio-padrão do efeito de lote
mu <- 0    # exp(mu)/(1+exp(mu)) = probabilidade de germinação
i <- 10    # níveis de lote
j <- 4     # número de repetições por lote
m <- 50    # número de sementes por ensaio
b <- rnorm(i, 0, sigma) # efeitos dos lotes
z <- gl(i, j)            # vetor de níveis de lote
y <- rbinom(i*j, size=m, prob=1/(1+exp(-mu-rep(b, e=j)))) # germinação observada

require(lattice)
xyplot(y~z, jitter.x=TRUE, col=z)

#-----------------------------------------------------------------------------
# função de log-verossimilhança

L <- function(par, y, z, m=50){
  ## par: vetor intercepto e parâmetro de variância (numeric)
  ## y: vetor valores observados para o número de sucessos (numeric)
  ## z: vetor de níveis de lote (factor)
  ## m: escalar número de sementes no ensaio (integer)
  ##-------------------------------------------------------
  ## lista em que cada slot são dos dados de um lote
  data <- split(y, z)
  ##-------------------------------------------------------
  ## verossimilhança para cada lote: f(y|b)*f(b)
  ## tem que função vetorial para usar na integrate()
  Li <- function(b, par, yi, m){
    sapply(b,
           function(bj){
             nu <- 1/(1+exp(-(par[1]+bj)))
             lj <- sum(dbinom(yi, prob=nu, size=m, log=TRUE))+
               dnorm(bj, mean=0, sd=exp(par[2]), log=TRUE)
             return(exp(lj))
           })
  }
  ##-------------------------------------------------------
  ## valor das i integrais
  int <-
    lapply(data,
           function(i){
             integrate(Li, -Inf, Inf, par=par,
                       yi=i, m=m)$value
           })
  ##-------------------------------------------------------
  ## valor da log-verossimilhança
  ll <- sum(log(unlist(int)))
  print(c(ll, par))
  return(ll)
  ##-------------------------------------------------------
}

#-----------------------------------------------------------------------------
# usa optim() para encontrar as estimativas de máxima verossimilhança

op <- optim(c(0,0), L, y=y, z=z, m=50,
            method="BFGS", control=list(fnscale=-1))
op$par   # estimativas de máxima verossimilhança para b0 e log(sigma)
op$value # valor do máximo da verossmilhança

#-----------------------------------------------------------------------------
# perfil de verossimilhança conjunto para b0 e log(sigma)

mu.grid <- seq(-2, 1, l=50)
lsigma.grid <- seq(-1, 1, l=50)
grid <- expand.grid(mu=mu.grid, lsigma=lsigma.grid)
grid$ll <- apply(grid, 1, L, y=y, z=z, m=m)

qch <- qchisq(c(0.9,0.95,0.99), df=2)
grid$dev <- -2*(grid$ll-op$value)

wireframe(dev~lsigma+mu, data=grid)

png("f022.png", w=500, h=500)
levelplot(dev~lsigma+mu, data=grid,
          xlab=expression(log(sigma)), ylab=expression(mu),
          panel=function(..., at, contour, region){
            panel.levelplot(..., at=at)#, at=94:195, contour=FALSE, region=TRUE)
            panel.contourplot(..., at=qch, contour=TRUE, region=FALSE)
            panel.abline(v=op$par[2], h=op$par[1], lty=3)
            panel.points(log(sigma), mu, pch=3, col=1)
          })
dev.off()

#-----------------------------------------------------------------------------
# usando a glmer()

require(lme4)

mm0 <- glmer(cbind(y, 50-y)~(1|z), family=binomial)
summary(mm0)
c(op$par[1], exp(op$par[2]))

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