########################################################################################## ## Exemplo de estimação em modelos Poisson ############################################### ########################################################################################## ## Modelo 1: ## Y_i ~ P(lambda) y <- rpois(100, lambda = 10) ## Funcao de verossimilhanca ll.pois <- function(lambda, y){ ll <- sum(dpois(y, lambda = lambda, log=TRUE)) return(ll) } ## Sequencia de lambda's lambda.seq <- seq(5, 15, l=100) ## Avaliando a funcao em cada ponto ll <- c() for(i in 1:100){ ll[i] <- ll.pois(lambda = lambda.seq[i], y=y) } ## outra forma evitando o for() lla <- sapply(lambda.seq, ll.pois, y=y) all.equal(ll, lla) rm(lla) ## Mostrando no grafico plot(ll ~ lambda.seq, type="l") ## Otimizando a funcao via optim estimativa <- optim(par = 10, fn = ll.pois, y =y, method="BFGS", control=list(fnscale=-1), hessian=TRUE) estimativa$par abline(v=estimativa$par) ## comparando com a solucao analitica disponível neste caso mean(y) ## variância da estimativa - obtida numéricamente ## hessiano numérico como aproximação da Informação (observada) (v.lambda <- (-estimativa$hessian)^-1) ## que **neste caso** é muito próxima à variância assintótica mean(y)/length(y) ## intervalo de confiança (assintótico/wald - aprox. quadratica) estimativa$par + c(-1, 1) * qnorm(0.975)*sqrt(v.lambda) ## Teste de razao de verossimlhanca -2*(ll.pois(8, y=y) - estimativa$value) qchisq(0.95, df=1) ########################################################################################## ## Modelo de regressao Poisson ########################################################### ########################################################################################## ## Modelo 2: ## Y_i ~ P(lambda_i) ## lambda_i = exp{beta_0 + beta_1 x} ## Criando a covariavel x <- seq(0,10, l=100) ## Fixando os valores dos parametros b0 = 1 b1 = 0.4 ## Calculando o eta (preditor linear) eta = b0 + b1*x ## Calculando o lambda lambda = exp(eta) ## Gerando a amostra y <- rpois(100, lambda=lambda) ## Plot contra a covariavel plot(y ~ x) plot(log(y) ~ x) ## Escrevendo a funcao de verossimilhanca formula <- "y ~ x" dados <- data.frame(y=y, x=x) ll <- function(par, Y, X){ eta <- X%*%par lambda = exp(eta) ll.pois <- sum(dpois(Y, lambda = lambda, log=TRUE)) return(ll.pois) } ll.reg.Poisson <- function(formula, data){ modelo <- model.frame(as.formula(formula), data=data) Y <- model.response(modelo) X <- model.matrix(as.formula(formula), data=data) inits <- rep(0, dim(X)[2]) estimativa = optim(par = inits, ll, Y = Y, X = X, control=list(fnscale =-1), method="BFGS",hessian=TRUE) pontual <- estimativa$par logLik <- estimativa$value std.error <- sqrt(diag(solve(-estimativa$hessian))) inf <- pontual - qnorm(0.975)*std.error sup <- pontual + qnorm(0.975)*std.error saida <- list(pontual = pontual, logLik = logLik, limite = cbind(inf, sup)) return(saida) } z <- seq(10, 0, l=100) data = data.frame(y, x, z) ll.reg.Poisson(y ~ x, data = dados) fit.glm <- glm(y ~ x, family=poisson) logLik(fit.glm) coef(fit.glm) confint.default(fit.glm) confint(fit.glm) ########################################################################################## ## Modelo Poisson com intercepto aleatório ############################################### ########################################################################################## ## Modelo 3: ## Y_i|b_i ~ P(lambda_i) ## lambda_i = exp{beta_0 + b_i} ## b_i ~ N(0, sigma^2) ## Fixando valores para os parametros b0 = 5 sigma <- 0.5 ## Amostrando os efeitos aleatorios bi <- rnorm(100, 0 , sd = sigma) ## calculando os lambda's lambda.i = exp(b0 + bi) ## Amostrando a resposta y = rpois(100, lambda = lambda.i) ### Escrevendo a funcao integrando integrando <- function(bi, b0, sigma, Y){ lambda <- exp(b0 + bi) int <- dpois(Y, lambda = lambda, log=TRUE) + dnorm(bi, mean = 0, sd = sigma, log=TRUE) return(exp(int)) } ## visualizando a função a ser integrada (para 1a observação) grid.b <- seq(-1,2, l = 100) ll.b <- c() for(j in 1:100){ ll.b[j] <- integrando(grid.b[j], b0 = 5, sigma = 0.5, Y = y[1]) } plot(ll.b ~ grid.b, type="l") ## na estimação, esta função vai ser: ## - integrada para cada observação ## - maximizada em relação aos parâmetros b0 e sigma ## visualizando novamente a função a ser integrada (para 10 observações) grid.b <- seq(-2,2, l = 100) ll.b <- matrix(0, 100, 10) for(i in 1:10){ for(j in 1:length(grid.b)){ ll.b[j,i] <- integrando(grid.b[j], b0 = 5, sigma = 0.5, Y = y[i]) } } matplot(grid.b, ll.b, type="l", lty=1) ## Resolvendo a integral via integrate ll.glmm <- function(par, Y){ b0 <- par[1] sigma <- par[2] ll.integral <- c() for(i in 1:length(Y)){ ll.integral[i] <- integrate(integrando, b0 = par[1], sigma=par[2], Y = y[i], lower=-Inf, upper=Inf)$value } ll.total <- sum(log(ll.integral)) return(ll.total) } ## Maximizando a funcao de log-verossimilhanca marginal estima.glmm <- optim(par=c(5, 0.5), ll.glmm, Y = y, control=list(fnscale=-1), method="L-BFGS-B", lower=c(-Inf, 1e-32), upper=c(Inf, Inf),hessian=TRUE) estima.glmm$par (std.error <- sqrt(diag(solve(-estima.glmm$hessian)))) estima.glmm$par + matrix(c(-1, -1, 1, 1), nc=2) * qnorm(0.975)*std.error ########################################################################################## ### Usando o pacote bbmle ################################################################ ########################################################################################## ll.bbmle <- function(b0, sigma, Y){ ll <- ll.glmm(par = c(b0, sigma), Y = Y) return(-ll) } ########################################################################################## ## Carregando o pacote bbmle ############################################################# ########################################################################################## require(bbmle) estimativa <- mle2(ll.bbmle, start=list(b0=5, sigma=0.5), method="L-BFGS-B", lower=list(b0 = -Inf, sigma = 1e-32), upper=list(b0 = Inf, sigma = Inf), data=list(Y=y)) estimativa ## Intervalos de confianca (intervalo <- confint(estimativa)) ## Perfil de verossimilhanca perf <- profile(estimativa) plot(perf)