########################################################################################## ## Modelo de regressão Gamma ############################################################# ########################################################################################## Y|x ~ G(mu, phi) log(mu) = beta0 + beta1*x ## Definindo as quantidades conhecidas beta0 <- 2 beta1 <- 3 phi = 10 x <- seq(0,1,l=100) mu <- exp(beta0 + beta1*x) ## Simulando do modelo set.seed(123) y <- rgamma(100, shape = phi, scale = mu/phi) ## Diagrama de dispersão plot(y ~x) ## Escrevendo a função densidade da gama na parametrização my.dgamma <- function(y, mu, phi, log=TRUE){ dens <- phi*log(phi) - lgamma(phi) - phi*log(mu) + (phi - 1)*log(y) - (phi*y)/mu if(log == FALSE){dens <- exp(dens)} return(dens) } ## Gráfico da função densidade de probabilidade curve(my.dgamma(x,mu = 10, phi = 10, log=FALSE),0,30) curve(my.dgamma(x,mu = 5, phi = 10, log=FALSE),0,30) curve(my.dgamma(x,mu = 10, phi = 20, log=FALSE),0,30) ## Função de ligação link <- function(mu){log(mu)} ## Inverso da função de ligação inv.link <- function(eta){exp(eta)} ## Função de log-verossimilhança logLik <- function(par, y, x){ beta0 <- par[1] beta1 <- par[2] phi <- par[3] eta <- beta0 + beta1*x mu <- inv.link(eta) ll <- sum(my.dgamma(y=y, mu = mu, phi = phi, log=TRUE)) return(ll) } ## Otimizando via optim fit1 <- optim(par=c(1,1,1), logLik, y=y, x=x, method="L-BFGS-B", lower=c(-Inf, -Inf, 1e-10), upper=c(Inf, Inf, Inf), control=list(fnscale=-1), hessian=TRUE) fit1$par ## Construindo intervalo Wald std.error <- sqrt(diag(solve(-fit1$hessian))) lim.sup <- fit1$par + qnorm(0.975)*std.error lim.inf <- fit1$par - qnorm(0.975)*std.error resul.optim <- cbind(fit1$par,std.error, lim.inf, lim.sup) resul.optim ## Visualizando o ajuste X <- model.matrix(~x) eta <- fit1$par[1] + fit1$par[2]*x mu <- inv.link(eta) plot(y ~ x) lines(x,mu) ## Avaliando o tempo computacional system.time(replicate(100,fit1 <- optim(par=c(1,1,1), logLik, y=y, x=x, method="L-BFGS-B", lower=c(-Inf, -Inf, 1e-10), upper=c(Inf, Inf, Inf), control=list(fnscale=-1), hessian=TRUE))) ########################################################################################## ## Obtendo a derivada analitica ########################################################## ########################################################################################## escore <- function(par, y, x){ beta0 <- par[1] beta1 <- par[2] phi <- par[3] eta <- beta0 + beta1*x mu <- inv.link(eta) dev.beta0 <- sum( ((phi*y)/mu) - phi) dev.beta1 <- sum(((phi*y*x)/mu) -phi*x) dev.phi <- sum(log(y) - y/mu + log(phi) - digamma(phi) - log(mu) + 1) esc <- c(dev.beta0,dev.beta1,dev.phi) return(esc) } escore(par = fit1$par, y=y,x=x) ## A mesma operação usando a derivada fit2 <- optim(par=c(1,1,1), logLik, gr = escore, y=y, x=x, method="L-BFGS-B", lower=c(-Inf, -Inf, 1e-10), upper=c(Inf, Inf, Inf), control=list(fnscale=-1), hessian=TRUE) fit2 std.error <- sqrt(diag(solve(-fit2$hessian))) lim.sup <- fit2$par + qnorm(0.975)*std.error lim.inf <- fit2$par - qnorm(0.975)*std.error resul2.optim <- cbind(fit2$par,std.error, lim.inf, lim.sup) resul2.optim resul.optim system.time(replicate(100,fit2 <- optim(par=c(1,1,1), logLik,gr=escore, y=y, x=x, method="L-BFGS-B", lower=c(-Inf, -Inf, 1e-10), upper=c(Inf, Inf, Inf), control=list(fnscale=-1), hessian=TRUE))) fit.glm <- glm(y ~ x, family=Gamma(link = "log")) coef(fit.glm) resul.optim system.time(replicate(100,fit3 <- glm(y ~ x, family=Gamma(link = "log"))))