Archive

Posts Tagged ‘optim’

Moda de uma amostra: método do histograma vs método kernel

Moda de uma amostra aleatória de uma distribuição contínua obtida pelo método do histograma e da densidade kernel.

Uma pergunta bem comum nos fóruns de estatística, geralmente associado à alunos de graduação, é de alguma função no R para obter a moda de um conjunto de dados. Bem, como moda se refere ao valor de maior frequência na amostra, a tarefa é simples para dados discretos: basta obter a distribuição de frequência dos valores (função table()) e tomar aquele associado à maior frequência.

Para dados contínuos, teoricamente, não teríamos dois valores observados iguais em uma amostra. Em amostras reais às vezes temos, seja pelo fato da amostra ser grande e/ou pelo instrumento de medida não ser preciso para diferenciar duas observações. Por exemplo, uma balança que registra massa em quilos vai estar sempre arredondando (censurando) a parte decimal dos valores. Então na prática até que podemos ter, eventualmente, valores repetidos em amostras de variáveis contínuas. Porém, não é assim que calculamos moda nessas amostras, porque valores repetidos podem não ocorrer e quando ocorrem, raramente excede à duas ocorrência de um mesmo valor.

A obtenção da moda pelo método do histograma é simples e está bem documentada nos livros texto da disciplina de estatística básica. Consiste em agrupar os dados em classes e fazer uma interpolação linear que considera a frequência da classe modal e de suas vizinhas. A interpolação linear usa semelhança de triângulos para encontrar a moda da amostra. A fórmula é essa

mo= l_i+\displaystyle A_c \left(\frac{f_m-f_e}{2f_m-f_e-f_d}\right )

em que l_i é o limite inferior da classe modal, A_c é a amplitude da classe modal, f_m é a frequência da classe modal, f_e é a frequência da classe à esquerda da modal, f_d é a frequência da classe à direita da modal.

Outra forma de obter a classe modal é calculando o valor da amostra que corresponde ao máximo da função de densidade kernel. É um método mais moderno e computacionalmente mais trabalhoso, pois deve-se ajustar a densidade kernel, fazer aproximação da função kernel por uma função interpoladora e depois obter o máximo dessa função. Mas como somos usuários de R, já temos há um bom tempo as ferramentas necessárias para isso.

No CMR abaixo eu implementei os métodos para obtenção da moda de uma amostra pelos dois métodos: do histograma e da densidade kernel. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição discreta

x <- rpois(100, lambda=12)
plot(sort(x))
tb <- table(x); tb
tb[which.max(tb)]

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição contínua

x <- rbeta(1000, 0.5,0.5)
ht <- hist(x) # histograma

#-----------------------------------------------------------------------------
# função para obter a moda pelo método do histograma, moda obtida por
# interpolação linear, baseado em semelhaça de triângulos

moda.hist <- function(ht, plotit=TRUE){
  ## ht: um objeto do uso da função hist()
  mcl <- which.max(ht$counts)
  li <- ht$breaks[mcl]
  width <- diff(ht$breaks[mcl+0:1])
  counts <- c(0,ht$counts,0)
  delta <- abs(diff(counts[1+mcl+(-1:1)]))
  moda <- li+width*delta[1]/sum(delta)
  cols <- rep(5, length(ht$counts))
  cols[mcl] <- 3
  if(plotit==TRUE){
    plot(ht, col=cols)
    abline(v=moda)
  }
  return(moda=moda)
}

ht <- hist(x)
moda.hist(ht)

#-----------------------------------------------------------------------------
# função para obter a moda da função densidade kernel

moda.dens <- function(dn, plotit=TRUE){
  ## dn: um objeto do uso da função density()
  ini <- dn$x[which.max(dn$y)]
  fx <- approxfun(dn$x, dn$y)
  op <- optim(c(ini), fx, method="BFGS", control=list(fnscale=-1))
  if(plotit==TRUE){
    plot(dn)
    abline(v=op$par, col=2)
  }
  return(moda=op$par)
}

dn <- density(x, kernel="triangular", bw=0.01)
moda.dens(dn)

dn <- density(x, kernel="gaussian", bw=0.05)
moda.dens(dn)

#-----------------------------------------------------------------------------
# gráfico

png("f027.png", 400, 300)
x <- rgamma(1000, 5, 5)
ht <- hist(x, freq=FALSE, col="gray86")
mo1 <- moda.hist(ht, plotit=FALSE)
abline(v=mo1, col=2, lwd=2)
dn <- density(x)
lines(dn, col=4, lwd=2)
mo2 <- moda.dens(dn, plotit=FALSE)
abline(v=mo2, col=3, lwd=2)
box()
dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, , ,

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

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