Arquivo

Posts Tagged ‘verossimilhança’

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!

Anúncios

Reparametrização do modelo quadrático para ponto crítico

Perfil de log-verossimilhaça representada pelo módulo da raíz da função deviance para cada um dos parâmetros do modelo quadrático reparametrizado. Linhas tracejadas indicam os limites dos intervalos de confiança de 99, 95, 90, 80 e 50%. Padrões de perfil diferentes de um V indicam intervalos assimétricos.

O emprego de funções polinomiais é muito frequente em ciências agrárias. O polinômio de segundo grau ou modelo quadrático é a escolha mais natural e imediata de modelo para quando a variável dependente apresentam uma tendência/relação não linear com a variável independente. Não raramente, o emprego desse modelo motiva também a estimação do ponto crítico (mínimo ou máximo). Para isso, o que se faz é estimar os parâmetros do modelo e via regras do cálculo obter a estimativa pontual do ponto crítico. O problema é que normalmente se necessita e não se faz, por desconhecimento de métodos, associação de incerteza à essa estimativa. É sobre isso que vou tratar agora: como associar incerteza ao ponto crítico de um modelo quadrático.

O modelo quadrático tem essa equação

m(x) = a+b\cdot x + c\cdot x^2,

em que a é o intercepto, ou seja, E(Y|x=0) = ab é valor da derivada do modelo na origem, ou seja, d(m_0(x))/dx = b, e c mede a curvatura. Derivando o modelo, igualando a zero e resolvendo, obtemos que o ponto crítico é

x_m = -b/(2\cdot c),

portanto, x_m é uma função das estimativas dos parâmetros, e como tal, deve apresentar incerteza associada.

Uma forma de associar incerteza é empregar o método delta para encontrar erros padrões. No R, o mesmo está disponível em msm::deltamethod(). Outra forma é usar um simulação boostratp e aplicar nas amostras bootstrap a relação acima. E o terceiro, e que eu mais gosto pois tem uma série de vantagens, é usar o perfil de log-verossimilhança.

Para usar o perfil de verossimilhança eu tive que escrever o modelo de forma que o ponto crítico fosse um parâmetro dele. Reparametrizar um modelo nem sempre é uma tarefa fácil. Nesse caso o procedimento escrito é curto e simples de entender. Pórem eu demorei cerca de uma hora com tentativa-erro até usar a seguinte abordagem: imagine uma função que tenha ponto de máximo. Podemos aproximar essa função ao redor de um valor x_0 por série de Taylor. Para uma aproximação de segundo grau temos

f(x) \approx f(x_0)+f'(x_0)(x-x_0)+0.5\cdot f''(x_0)(x-x_0)^2.

Agora imagine que x_0 é um ponto crítico da função. O termo linear da aproximação é 0 no ponto crítico por definição. Se a função que estamos aproximando for perfeitamente quadrática, esse procedimento nos leva a reparametrização que desejamos. Assim, o modelo quadrático reparametrizado é

n(x) = y_m+c\cdot (x-x_m)^2.

em que x_m é o nosso parâmetro de interesse, o ponto crítico, y_m é o valor da função no ponto crítico, ou seja, E(Y|x=x_m) = y_m, c é a métade da segunda derivada da função no ponto crítico, ou seja, c=0.5\cdot f''(x_m), é uma curvatura. Então nos reparametrizamos o model quadrático para ter o ponto crítico como parâmetro. Nessa reparametrização obtvemos 3 parâmetros com interpretação simples. Em termos de interpretação, se c=0 não temos curvatura, portanto, não temos ponto crítico e devemos ajustar um modelo linear da reta. Se temos curvatura podemos estimar o ponto crítico e o valor da função no ponto crítico. Então imagine que você estuda doses de nutrientes para uma cultura, não é fundamental saber qual a dose que confere o máximo e qual a produção esperada nessa dose? Melhor que isso é poder associar incerteza à essas quantidades.

O modelo reparametrizado, diferente do original, é um modelo não linear. Por isso, para estimar os parâmetros no R usaremos a função nls() e não a lm(). (Visite os posts passados para ver mais coisas relacionadas à modelos não lineares.) Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile(). Consulte a documentação dessas funções para saber o que elas de fato fazem.

No CMR abaixo eu ilustro a aplicação desses métodos. É recomendando ao leitor que compare os intervalos de confiança, que repita os procedimento simulando os dados com outros parâmetros, tamanhos e espaços, e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos (teste de hipótese) e o modelo quadrático para duas regressoras (bidimensional). Até a próxima rídicula.

#-----------------------------------------------------------------------------
# dados simulados

set.seed(23456)
x <- 1:12
y <- 3+0.17*x-0.01*x^2+rnorm(x,0,0.05)
plot(x, y)                         # dados observados
curve(3+0.17*x-0.01*x^2, add=TRUE) # curva que gerou os dados

#-----------------------------------------------------------------------------
# ajuste do y = a+b*x+c*x², é um modelo linear

m0 <- lm(y~x+I(x^2))
summary(m0)

#-----------------------------------------------------------------------------
# método delta

require(msm)
coef <- coef(m0)[2:3]; vcov <- vcov(m0)[2:3,2:3]
xm <- -0.5*coef[1]/coef[2]; xm                         # estimativa
sd.xm <- deltamethod(~(-0.5*x1/x2), coef, vcov); sd.xm # erro padrão
xm+c(-1,1)*qt(0.975, df=df.residual(m0))*sd.xm         # IC 95%

#-----------------------------------------------------------------------------
# simulação (bootstrap paramétrico)

require(MASS)
aa <- mvrnorm(1000, mu=coef, Sigma=vcov)
xm.boot <- apply(aa, 1, function(i) -0.5*i[1]/i[2])
mean(xm.boot)                      # média bootstrap
sd(xm.boot)                        # desvio-padrão bootstrap
quantile(xm.boot, c(0.025, 0.975)) # IC boostrap

#-----------------------------------------------------------------------------
# ajuste do y = ym+c*(x-xm)², é um modelo não linear

plot(x, y)
abline(h=3.7, v=8.5, lty=3) # aproximação visual do chute
start <- list(ym=3.7, xm=8.5, c=coef(m0)[3])
max(fitted(m0))             # boa aproximação de chute
x[which.max(fitted(m0))]    # boa aproximação de chute
n0 <- nls(y~ym+c*(x-xm)^2, start=start)
summary(n0) # estimativa do xm que é parâmetro do modelo, estimado diretamente

#-----------------------------------------------------------------------------
# gráfico dos observados com a curva ajustada

plot(x, y)
with(as.list(coef(n0)),{
  curve(ym+c*(x-xm)^2, add=TRUE)
  abline(h=ym, v=xm, lty=3)
})

#-----------------------------------------------------------------------------
# veja a igualdade dos modelos

cbind(coef(m0),
      with(as.list(coef(n0)), c(ym+c*xm^2, -2*c*xm, c)))
deviance(m0)
deviance(n0)
plot(fitted(m0), fitted(n0)); abline(0,1)

#-----------------------------------------------------------------------------
# intervalos de confiança

confint.default(n0) # IC assintóticos, são simétricos por construção
confint(n0)         # IC baseado em perfil de log-verossimilhança

#-----------------------------------------------------------------------------
# gráficos de perfil de verossimilhança

#png(file="f030.png", width=500, height=175)
par(mfrow=c(1,3))
plot(profile(n0))
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# representação do ajuste

pred <- data.frame(x=seq(0,13,l=30))
pred <- cbind(as.data.frame(predict(m0, newdata=pred, interval="confidence")),
              pred)
str(pred)
ic <- confint(n0); ic

plot(x, y, type="n")              # janela vazia
with(pred,                        # bandas de confiança
     polygon(c(x,rev(x)), c(lwr,rev(upr)), col="gray90", border=NA))
points(x, y)                      # valores observados
with(pred, lines(x, fit))         # valores preditor
abline(v=ic[2,], h=ic[1,], lty=2) # IC para xm e ym

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

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]))

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