|
Métodos Computacionais para Inferência Estatística |
Modelo Poisson com efeito aleatório
Prof. Wagner Hugo Bonat LEG/UFPR
2019
1 Introdução
A distribuição de Poisson é uma escolha comum para a modelagem de dados de contagem. No entanto é também bem documentado na literatura que tal escolha é muito restritiva em muitas situações reais, devido a forte suposição de equidispersão, ou seja, esperança igual a variância inerente a distribuição de Poisson.
Sendo \(Y\) uma v.a com distribuição Poisson sua função de probabilidade é dada por
\[P(Y = y) = \exp\{-\lambda\} \frac{\lambda^y}{y!},\] onde \(\lambda > 0\) e \(y \in N\). Nesta notação temos que \(\mathrm{E}(Y) = \lambda\) e \(\mathrm{Var}(Y) = \lambda\) evidenciando a restrição de igualdade entre esperança e variância. Uma forma de lidar com superdispersão, ou seja, \(\mathrm{Var}(Y) > \mathrm{E}(Y)\) é incluir uma nova fonte de variação na distribuição de Poisson. Essa nova fonte de variação combinada com a distribição de Poisson irá gerar uma nova distribuição em que potencialmente a \(\mathrm{V}(Y) > \mathrm{E}(Y)\) e portanto será adequada para lidar com superdispersão.
Uma forma de induzir tal acréscimo de variabilidade é introduzir heterogeineidade entre as v.a, ou seja, ao invés de todas as \(Y_i\) v.a’s terem o mesmo parâmetro \(\lambda\), nos podemos assumir que cada v.a \(Y_i\) tem distribuição Poisson com um parâmetro \(\lambda_i\) específico. O problema agora é como especificar um parâmetro para cada v.a e obviamente estimar o modelo.
A idéia é decompor o parâmetro \(\lambda\) da distribuição Poisson em dois componentes, um comum a todas as v.a e outro que varia v.a por v.a. E para esta variação extra nos dizemos que ela também vem de uma outra distribuição de probabilidade. Uma forma comum de fazer é assumindo a distribuição normal. Neste caso \(g(\lambda_i) = \alpha + b_i\), onde \(b_i \sim N(0, \sigma^2)\) e \(g:\Re_{+} \to \Re\). Note que sendo \(b_i\) normal seu suporte é os reais, porém \(\lambda_i\) deve ser estritamente positivo, assim precisamos especificar \(g^{-1}(\cdot)\) de tal forma que receba um real e retorne um real positivo. Uma função muito popular com esta propriedade é a função exponencial. A ideia é exatamente a mesma que uma função de ligação no contexto de modelos lineares generalizados.
Outro ponto importante é que os desvios \(b_i\) não são observados baseado nos dados. Eles podem ser interpretados com condições que levam a heterogeneidade não observada entre as observações que causam uma variabilidade maior do que a esperada sob o modelo Poisson. Em geral chama-se os \(b_i's\) de efeitos aleatórios ou também de efeitos latentes.
Usando estas suposições nosso modelo fica especificado como \[Y_i | b_i \sim P(\lambda_i), \\ b_i \sim N(0, \sigma^2) \quad \text{onde} \quad \lambda_i = \exp(\alpha + b_i). \] Note que nosso modelo agora tem dois parâmetros \(\boldsymbol{\theta} = (\alpha, \sigma^2)^{\top}.\)
Nosso objetivo é estimar o vetor de parâmetro e fornecer inferência, em geral através de intervalos de confiança e testes de hipóteses. Para isto podemos adotar o método de máxima verossimilhança.
2 Função de verossimilhança
Conforme discutido durante o curso a função de verossimilhança é definida como a distribuição conjunta de todas as variáveis aleatórias envolvidas no modelo, porém vista como função dos parâmetros. Neste exemplo, temos duas variáveis aleatórias \(Y_i|b_i\) e \(b_i\), porém \(b_i\) não é observado. Neste caso precisamos obter a distribuição marginal de \(Y_i\). Usando resultados básicos de probabilidades, sabemos que a distribuição marginal de \(Y_i\) é obtida integrando \(b_i\) na distribuição conjunta de \((Y_i,b_i)\), ou seja, \[f(y_i|\boldsymbol{\theta}) = \int f(y_i|b_i) f(b_i) db_i.\] Sendo assim, a função de verossimilhança é dada por \[L(\boldsymbol{\theta}|y_i) = \int f(y_i|b_i) f(b_i) db_i.\] Assumindo independência entre as observações dado os efeitos aleatórios, podemos escrever a verossimilhança para o nosso modelo como sendo \[ L(\boldsymbol{\theta}|y_i) = \prod_{i=1}^n \int_{-\infty}^{\infty} \exp\{-\lambda_i\} \frac{\lambda_i^{y_i}}{y_i!} \frac{1}{\sqrt{2\pi}\sigma} \exp\left \{ \frac{1}{2\sigma^2} b_i^2 \right \}. \]
Note que para avaliar a função de verossimilhança precisamos resolver \(n\) integrais unidimensionais, onde o integrando é o produto de uma distribuição de Poisson por uma normal padrão. Obviamente precisamos de métodos numéricos para resolver tal integral.
2.1 Avaliando a função de verossimilhança
Nosso primeiro objetivo é resolver numéricamente a integral dentro da função de verossimilhança. Note que o dominio de integração é os números reais, assim temos pelo menos três diferentes formas de resolver tal integral: Gauss-Hermite, Laplace e Monte Carlo. Nos slides da aula foram apresentadas funções genéricas que implementam cada um destes métodos as quais reproduzimos abaixo:
require(pracma)
# Loading required package: pracma
# Gauss-Hermite
gauss_hermite <- function(integrando, n.pontos, ...) {
pontos <- gaussHermite(n.pontos)
integral <- sum(pontos$w*integrando(pontos$x,...)
/exp(-pontos$x^2))
return(integral)
}
# Laplace
laplace <- function(funcao, otimizador,n.dim, ...) {
integral <- -999999
inicial <- rep(0,n.dim)
temp <- try(optim(inicial, funcao,..., method=otimizador,
hessian=TRUE, control=list(fnscale=-1)))
if(class(temp) != "try-error"){
integral <- exp(temp$value) * (exp((n.dim/2)*log(2*pi) -
0.5*determinant(-temp$hessian)$modulus))
}
return(integral)
}
# Monte-Carlo
monte.carlo <- function(funcao, n.pontos, ...) {
pontos <- rnorm(n.pontos)
norma <- dnorm(pontos)
integral <- mean(funcao(pontos,...)/norma)
return(integral)
}Em todos os casos precisamos implementar a função integrando. Por questões de precisão computacional recomenda-se implementar a função sempre em escala log e depois retornar a escala de interesse.
integrando <- function(b, theta, y, log = FALSE) {
out <- dpois(x = y, lambda = exp(theta[1] + b), log = TRUE) +
dnorm(b, mean = 0, sd = sqrt(theta[2]), log = TRUE)
if(log == FALSE) {out <- exp(out)}
return(out)
}Vamos fazer um gráfico da função que queremos integrar.
b = seq(-3,3, l = 100)
fb <- integrando(b = b, theta = c(log(10), 0.25), y = 10)
plot(fb ~ b, type = "l")
Podemos tentar resolver a integral numericamente para um dado vetor de parâmetros e um valor observado de \(y\). Note que para a integração via aproximação de Laplace o integrando tem que ser passado em escala log. Para conferir nossas implementações, vamos também resolver a integral usando a função integrate() do R.
theta = c(log(10), 0.25)
y = 10
# Gauss-Hermite
gauss_hermite(integrando = integrando, n.pontos = 100,
theta = theta, y = y)
# [1] 0.06695328
# Laplace
laplace(funcao = integrando, otimizador = "BFGS", n.dim = 1,
theta = theta, y = y, log = TRUE)
# [1] 0.06687412
# attr(,"logarithm")
# [1] TRUE
# Monte Carlo
monte.carlo(funcao = integrando, n.pontos = 10000, theta = theta,
y = y, log = FALSE)
# [1] 0.0668426
# Nativa do R
integrate(f = integrando, lower = -Inf, upper = Inf, theta = theta, y = y)
# 0.06695328 with absolute error < 4.4e-062.2 Número de pontos de integração
Uma pergunta natural em termos práticos é quantos pontos de integração utilizar nos algoritmos de Gauss-Hermite e Monte Carlo. Essa pergunta é extremamente relevante, porém não existe uma resposta que funcione em todos os casos. O número de pontos necessários vai depender do modelo e também dos valores dos parâmetros e das observações. Uma estratégia comum é avaliar a integral para um número crescente de pontos de integração e verificar quando o valor estabiliza. Veja o código abaixo
n = seq(10, 200, by = 5)
vv <- c()
for(i in 1:length(n)) {
vv[i] <- gauss_hermite(integrando = integrando, n.pontos = n[i],
theta = theta, y = y)
}
plot(vv ~ n, type = "o")
Pelo gráfico parece claro que após 20 pontos de integração o valor da integral estabilizou. Então este seria um valor razoável, ao menos para esta configuração de parâmetros e valor observado.
Um procedimento similar pode ser usado para definir o número de simulações Monte Carlo necessário para aproximar a integral.
n = seq(10, 20000, by = 100)
vv <- c()
for(i in 1:length(n)) {
vv[i] <- monte.carlo(funcao = integrando, n.pontos = n[i],
theta = theta, y = y)
}
plot(vv ~ n, type = "o")
Fica claro na Figura acima que um número muito maior de avaliações da função é necessário no método de Monte Carlo. Além disso, mesmo com 20000 simulações o valor da integral ainda sofre pequenas variações.
Finalmente podemos escrever uma função para avaliar a distribuição marginal de \(Y\) usando algum dos nossos algoritmos numéricos.
dpois_gauss <- function(y, theta, method, n_points = NULL) {
if(method == "GH") {
out <- gauss_hermite(integrando = integrando, n.pontos = n_points,
theta = theta, y = y)
}
if(method == "LAPLACE") {
out <- laplace(funcao = integrando, otimizador = "BFGS", n.dim = 1,
theta = theta, y = y, log = TRUE)
}
if(method == "Monte Carlo") {
out <- monte.carlo(funcao = integrando, n.pontos = n_points,
theta = theta, y = y)
}
if(method == "integrate") {
out <- integrate(f = integrando, lower = -Inf, upper = Inf,
theta = theta, y = y)
}
return(out)
} Podemos fazer uma gráfico desta distribuição ma,rginal para alguns valores de \(\alpha\) e \(\sigma^2\).
y_grid <- 0:100
dd1 <- c()
dd2 <- c()
dd3 <- c()
dd4 <- c()
for(i in 1:length(y_grid)) {
dd1[i] <- dpois_gauss(y = y_grid[i], theta = c(log(10), 0.05),
method = "GH", n_points = 100)
dd2[i] <- dpois_gauss(y = y_grid[i], theta = c(log(10), 0.25),
method = "GH", n_points = 100)
dd3[i] <- dpois_gauss(y = y_grid[i], theta = c(log(10), 0.5),
method = "GH", n_points = 100)
dd4[i] <- dpois_gauss(y = y_grid[i], theta = c(log(10), 1),
method = "GH", n_points = 100)
}
par(mfrow = c(2,2), mar=c(2.6, 2.8, 1.2, 0.5), mgp = c(1.6, 0.6, 0))
plot(dd1 ~ y_grid, ylab = expression(P(Y == y)), xlab = "y",
type = "h", main = expression(sigma^2 == 0.05))
plot(dd2 ~ y_grid, ylab = expression(P(Y == y)), xlab = "y",
type = "h", main = expression(sigma^2 == 0.25))
plot(dd3 ~ y_grid, ylab = expression(P(Y == y)), xlab = "y",
type = "h", main = expression(sigma^2 == 0.50))
plot(dd4 ~ y_grid, ylab = expression(P(Y == y)), xlab = "y",
type = "h", main = expression(sigma^2 == 1))
Conforme gostaríamos o novo parâmetro \(\sigma^2\) controla a variabilidade da nova distribuição, e portanto temos um modelo em que a variância é maior do que a média.
3 Simulando do modelo
Uma vez que especificamos o modelo completamente a simulação de observações do modelo é trival, conforme mostra o código abaixo.
set.seed(123)
n = 100
b <- rnorm(n, mean = 0, sd = 0.25^2)
lambda <- exp(log(10) + b)
y <- rpois(n, lambda = lambda)
plot(table(y), ylab = expression(P(Y == y)), xlab = "y")
4 Estimação e Inferência
Dado um conjunto de \(n\) observações do modelo, podemos estimar os parâmetros \(\alpha\) e \(\sigma^2\) pelo método de máxima verossimilhança. Assumindo observações independentes a função de log-verossimilhança é simplesmente a soma do log das distribuições marginais vista como uma função dos parâmetros, ou seja,
ll <- function(theta, y, method = "GH", n_points = 100) {
out <- c()
for(i in 1:length(y)) {
out[i] <- log(dpois_gauss(y = y[i], theta = theta,
method = method, n_points = n_points))
}
return(sum(out))
}Avaliando a verossimilhança em alguns pontos para verificar a implementação.
ll(theta = c(log(10), 0.25), y = y)
# [1] -285.1257
ll(theta = c(log(10), 0.05), y = y)
# [1] -268.7464
ll(theta = c(log(10), 0.0001), y = y)
# [1] -6119.546O próximo passo consiste em maximizar a função de log-verossimilhança. Para isto vamos usar o algoritmo de Newton. Em sala foi fornecido o algoritmo Newton, aqui faremos algumas adaptações para que possamos usar os parâmetros extra da função de log-verossimilhança que controlam o algoritmo que resolve a integral dentro da verossimilhança.
newton <- function(fx, jacobian, x1, tol = 1e-04, max_iter = 10, ...) {
solucao <- matrix(NA, ncol = length(x1), nrow = max_iter)
solucao[1,] <- x1
for(i in 1:max_iter) {
J <- jacobian(solucao[i,], ...)
#J <- IO(solucao[i,], y = y)
grad <- as.numeric(fx(solucao[i,], ...))
#grad <- as.numeric(fx(solucao[i,], y = y))
solucao[i+1,] = solucao[i,] - solve(J, grad)
print(solucao[i+1,])
if( sum(abs(solucao[i+1,] - solucao[i,])) < tol) break
}
return(solucao)
}O algoritmo de Newton exige o cálculo da primeira e segunda derivada da função objetivo. Note que neste caso tais derivadas são impossíveis de serem obtidas, uma vez que a função de log-verossimilhança é dada por uma integral que não tem solução analitica. Sendo assim, vamos usar os métodos numéricos que vimos em sala para obter uma aproximação para a primeira e a segunda derivada da função de log-verossimilhança. Para facilitar vamos usar as funções do pacote rootSolve.
## Função escore
fx <- function(theta, y, method = "GH", n_points = 100) {
out <- gradient(f = ll, x = theta, y = y, method = method, n_points = n_points)
return(out)
}
## Negativo da Informação observada
IO <- function(theta, y, method = "GH", n_points = 100) {
out <- pracma::hessian(f = ll, x0 = theta, y = y,
method = method, n_points = n_points)
return(out)
}Resolvendo o sistema não-linear usando o método de Newton
require(rootSolve)
require(pracma)
fit <- newton(fx = fx, jacobian = IO, x1 = c(log(10), 0.25^2),
y = y, max_iter = 10)
# [1] 2.248236089 0.009089622
# [1] 2.26898613 0.01149581
# [1] 2.26797562 0.01410513
# [1] 2.26657829 0.01708982
# [1] 2.26460873 0.02097707
# [1] 2.26159655 0.02688336
# [1] 2.25861808 0.03298173
# [1] 2.25802022 0.03436426
# [1] 2.25800670 0.03440058
# Solução
fit[dim(na.exclude(fit))[1],]
# [1] 2.25800670 0.03440058Usando a distribuição assintótica do MLE, podemos facilmente calcular intervalos de confiança.
Est <- fit[dim(na.exclude(fit))[1],]
I0_obs <- IO(theta = Est, y = y)
VarCov <- solve(-I0_obs)
Var <- sqrt(diag(VarCov))
data.frame("Estimates" = Est, "Std.Error" = Var,
"Z-value" = Est/Var)
# Estimates Std.Error Z.value
# 1 2.25800670 0.03817094 59.155122
# 2 0.03440058 0.01905075 1.805734| Métodos Computacionais para Inferência Estatística | Prof. Wagner Hugo Bonat |
