1. Simular datos para un modelo linear mixto con intercepto aleatorio

Se simula el modelo \(Y=\beta_0 + \beta_1X + b_0 + \epsilon\) donde \(b_0 \sim N(0, \sigma_{b_0}^2)\) y \(\epsilon \sim N(0, \sigma^2)\) para mgrupos con nindividuos cada uno.

interceptoNorm1 <-  function (mgrupos, nindividuos, beta0, beta1, sigma2, sigma2b0){

  grupo.Total <- factor(rep(1:mgrupos, each = nindividuos))
  n <- length(grupo.Total) #( es lo mismo que n = mgrupos*nindividuos )
  
  intcept.Aleatorio <- rnorm(mgrupos, sd = sqrt(sigma2b0)) # Este es bo intercepto aleatorio
  er <- rnorm(n, mean = 0, sd = sqrt(sigma2))  #Este es el vector de errores
  x <- runif(n, 0, 10)
  y <- beta0 + beta1 * x + intcept.Aleatorio[grupo.Total] + er
  aleatorio <- data.frame(aleatorio=intcept.Aleatorio[grupo.Total])
  
  res <- data.frame(x, y, grupo.Total,aleatorio)
  names(res) <- c("x", "y", "grupo.Total", "aleatorio")
 
   
  
 d <- data.frame(x, y, grupo.Total)
  library(ggplot2)
 a <- ggplot(d, aes(x, y, color = grupo.Total) ) + geom_point() + labs(colour = "Grupo / Cluster") + theme(plot.margin = unit(c(0.8, 0.5, 0.8, 0.5), "cm")) 
 
 print(a)
  return(res)
}

Especificamente se simula el modelo \(Y=2 + 4X + b_0 + \epsilon\) donde \(b_0 \sim N(0, 80)\) y \(\epsilon \sim N(0, 0.8)\) para 5 grupos con 800 individuos cada uno.

datos1 <- interceptoNorm1(5,800,2,4,0.8,80)

head(datos1)
##          x        y grupo.Total aleatorio
## 1 6.673652 41.08731           1  13.67703
## 2 7.250569 44.74809           1  13.67703
## 3 5.262907 36.93826           1  13.67703
## 4 1.045350 20.78237           1  13.67703
## 5 2.431469 25.71448           1  13.67703
## 6 8.487385 48.40468           1  13.67703

2. Estimar los parámetros (\(b_0\), \(b_1\), \(\sigma^2\), \(\sigma_{b_0}^2\)) con tres paquetes en R:

a) lme

library(nlme)
lmeprueba <- lme(y ~ x,
                      random = ~ 1|aleatorio,
                      data = datos1)
 
summary(lmeprueba)
## Linear mixed-effects model fit by REML
##  Data: datos1 
##       AIC      BIC    logLik
##   10581.5 10606.68 -5286.752
## 
## Random effects:
##  Formula: ~1 | aleatorio
##         (Intercept)  Residual
## StdDev:    8.780478 0.9007279
## 
## Fixed effects: y ~ x 
##                Value Std.Error   DF  t-value p-value
## (Intercept) 3.329903  3.926854 3994   0.8480  0.3965
## x           3.996196  0.004917 3994 812.6491  0.0000
##  Correlation: 
##   (Intr)
## x -0.006
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.39699430 -0.69820950 -0.02441304  0.67827422  3.51444594 
## 
## Number of Observations: 4000
## Number of Groups: 5
fixef(lmeprueba) # b0 y b1
## (Intercept)           x 
##    3.329903    3.996196
sigma(lmeprueba)^2 # sigma2
## [1] 0.8113108
getVarCov(lmeprueba, type = "random") #sigma2b0
## Random effects variance covariance matrix
##             (Intercept)
## (Intercept)      77.097
##   Standard Deviations: 8.7805
#methods(class = "lme")

b) lmer

library(lme4)
lmerprueba <- lmer( y ~ x + (1|aleatorio), data = datos1)
summary(lmerprueba)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | aleatorio)
##    Data: datos1
## 
## REML criterion at convergence: 10573.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3970 -0.6982 -0.0244  0.6783  3.5144 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  aleatorio (Intercept) 77.0969  8.7805  
##  Residual               0.8113  0.9007  
## Number of obs: 4000, groups:  aleatorio, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 3.329903   3.926856   0.848
## x           3.996196   0.004917 812.649
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.006
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
fixef(lmerprueba) # b0 y b1
## (Intercept)           x 
##    3.329903    3.996196
sigma(lmerprueba)^2 #sigma2
## [1] 0.8113108
as.numeric(VarCorr(lmerprueba)$aleatorio) # sigma2b0
## [1] 77.09688

c) glmer

library("lme4")
glmerprueba <- glmer( y ~ x + (1|aleatorio), family = gaussian, data = datos1)
summary(glmerprueba)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | aleatorio)
##    Data: datos1
## 
## REML criterion at convergence: 10573.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3970 -0.6982 -0.0244  0.6783  3.5144 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  aleatorio (Intercept) 77.0969  8.7805  
##  Residual               0.8113  0.9007  
## Number of obs: 4000, groups:  aleatorio, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 3.329903   3.926856   0.848
## x           3.996196   0.004917 812.649
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.006
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
fixef(glmerprueba)  # b0 y b1
## (Intercept)           x 
##    3.329903    3.996196
sigma(glmerprueba)^2 # sigma2
## [1] 0.8113108
as.numeric(VarCorr(glmerprueba)$aleatorio) # sigma2b0
## [1] 77.09688

3. Resumen

A continuación una tabla comparativa con los resultados obtenidos con las 3 funciones.

todos <- data.frame(c("lme", "lmer", "glmer", "valor real"), c(effijolme[1],effijolmer[1],effijoglmer[1],2), c(effijolme[2],effijolmer[2],effijoglmer[2],4), c(sigmalme,sigmalmer,sigmaglmer,0.8), c(sigmab0lme,sigmab0lmer,sigmab0glmer,80))

library(knitr)
kable(todos, col.names = c("Función" ,"$b_0$", "$b_1$", "$\\sigma^2$", "$\\sigma_{b_0}^2$"))
Función \(b_0\) \(b_1\) \(\sigma^2\) \(\sigma_{b_0}^2\)
lme 3.329903 3.996196 0.8113108 77.09680
lmer 3.329903 3.996196 0.8113108 77.09688
glmer 3.329903 3.996196 0.8113108 77.09688
valor real 2.000000 4.000000 0.8000000 80.00000