Más sobre MLE

Objetivos:

Familiarizarse con la formulación de funciones de likelihood
Simular datos con distintas funciones de likelihood
Ajustar los datos simulados usando mle2

Empecemos con una regresión lineal normal con una sola variable predictora. En este caso el modelo es:

$$ \begin{aligned} & y_i \sim N(\mu_i, \sigma^{2}) \\ & \mu_i = \beta_0 + \beta_1 \times x_i\\ & \textrm{para } i \textrm{ de } 1:n \\ \end{aligned} $$

Para simular datos de acuerdo con este modelo, primero tenemos que definir los parámetros y la variable predictora (también llamada covariable) y el número de observaciones

b0 <- 0
b1 <- 1
s <- 3
n <- 30
set.seed(1234)
x <- runif(n, 0, 30)

Ahora podemos calcular μ para cada x y simular datos con una distribución Normal

m <- b0 + b1 * x
y <- rnorm(n, mean = m, sd = s)

Este modelo se puede ajustar fácilmente en R usando el comando lm, pero vamos a tomarnos el trabajo de escribir la función de likelihood y ajustarlo con mle2

neg.log.lik.lm <- function(pars) {
    b0 <- pars[1]
    b1 <- pars[2]
    s <- pars[3]
    m <- b0 + b1 * x
    nll <- -sum(dnorm(y, mean = m, sd = s, log = TRUE))
    return(nll)
}

library(bbmle)

parnames(neg.log.lik.lm) <- c("b0", "b1", "s")

fit.lm <- mle2(neg.log.lik.lm, start = c(b0 = 0, b1 = 0, s = 1), data = list(y = y, 
    x = x))

summary(fit.lm)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = neg.log.lik.lm, start = c(b0 = 0, b1 = 0, s = 1), 
##     data = list(y = y, x = x))
## 
## Coefficients:
##     Estimate Std. Error z value     Pr(z)    
## b0 -3.421537   0.843266 -4.0575 4.960e-05 ***
## b1  1.154153   0.052676 21.9105 < 2.2e-16 ***
## s   2.434711   0.314320  7.7460 9.486e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 138.526

Podemos graficar los datos y el modelo ajustado. Y ya que estamos comparamos con el resultado de lm y ver que da lo mismo... También podemos agregar al gráfico la relación verdadera entre x y μ

plot(x, y)

xvec <- seq(min(x), max(x), length = 100)
lines(xvec, fit.lm@coef[1] + fit.lm@coef[2] * xvec, lwd = 3)
lines(xvec, b0 + b1 * xvec, lwd = 2, col = "gray")
abline(lm(y ~ x), lty = 2, col = 2, lwd = 3)

Una buena práctica cuando hacemos regresiones de este tipo es centrar y estandarizar las variables predictoras (discutir por qué!):

xs <- scale(x)
fit.lms <- mle2(neg.log.lik.lm, start = c(b0 = 10, b1 = 10, s = 1), data = list(y = y, 
    x = xs))

Ejercicios

1- Hacer simulaciones y análisis como los de arriba pero reemplazando la Normal por: + Poisson + Binomial Negativa

2- Reemplazar la función lineal por una de saturación y usar una Poisson para la parte estocástica del modelo.

sessionInfo()

## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.5 (Yosemite)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] bbmle_1.0.17 knitr_1.11  
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-33      digest_0.6.8         mime_0.4            
##  [4] MASS_7.3-44          grid_3.2.2           knitrBootstrap_1.0.0
##  [7] formatR_1.2          magrittr_1.5         evaluate_0.7.2      
## [10] stringi_0.5-5        rmarkdown_0.8        tools_3.2.2         
## [13] stringr_1.0.0        markdown_0.7.7       numDeriv_2014.2-1   
## [16] yaml_2.1.13          htmltools_0.2.6

Juan Manuel Morales. 6 de Septiembre del 2015. Última actualización: 2015-09-28