Problemas con el modelo mixto de simulación

Question

Así que era tener una buena noche tratando de simular algunos modelos jerárquicos por mí mismo, y la estimación de sus parámetros. Primer modelo traté de simular y de una estimación de $$y_{ij}=\beta_0+\varepsilon_{ij},$$ donde $b_0=\gamma_{00}+u_{0j}$. Así que en realidad lo que yo era la estimación de se $y_{ij}=\gamma_{00}+u_{0j}+\varepsilon_{ij}$. Me han generado algunos datos y parámetros estimados $\gamma_{00}$ y las variaciones de los términos de error y que todo estaba bien:

set.seed(1)
N <- 100
nj <- 100
g00 <- 10
e <- rnorm(N*nj)
j <- c(sapply(1:N, function(x) rep(x, nj)))
uj <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))
d <- data.frame(j, y=g00+uj+e)
library(nlme)
lme(y~1, data=d, random=~1|j)

Linear mixed-effects model fit by REML
  Data: d 
  Log-restricted-likelihood: -14520.94
  Fixed: y ~ 1 
(Intercept) 
   10.00215 

Random effects:
 Formula: ~1 | j
        (Intercept) Residual
StdDev:   0.7752422 1.012683

Number of Observations: 10000
Number of Groups: 100 

Entonces traté de modelos diferentes: $$y_{ij}=\beta_0+\beta_1 x_{ij}+\varepsilon_{ij},$$ donde$\beta_0=\gamma_{00}+u_{0j}$$\beta_1=\gamma_{10}+u_{1j}$. Así que he tenido que estimar la ecuación $$y_{ij}=\gamma_{00}+u_{0j}+(\gamma_{10}+u_{1j}) \cdot x_{ij}+\varepsilon_{ij}=\gamma_{00}+\gamma_{10}x_ij+u_{0j}+u_{1j}x_{ij}+\varepsilon_{ij}.$$ Yo hice lo mismo que había antes, pero esta vez lme no convergen. He intentado esto, pero no parece funcionar.

g10 <- 10
u0j <- uj
u1j <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))
x1 <- rnorm(N*nj)
d1 <- data.frame(j, y=g00+u0j+(g10+u1j)*x1, x1)
lme(y~1+x1, data=d1, random=~1+x1|j)

Aquí es lo que la última llamada a lme escupir:

Error in lme.formula(y ~ 1 + x1, data = d1, random = ~1 + x1 | j) : 
  nlminb problem, convergence error code = 1
  message = false convergence (8)

¿Qué me puede sugerir a hacer? Quizás el problema es que con mi modelo de especificación, redundante parámetros, singular de la matriz o de otra cosa?

Answer 1

Lamentablemente no puedo explicar donde está el problema, pero esto funciona:

lme(y~1+x1, data=d1, random=~1+x1|j, control = lmeControl(opt = "optim"))
Linear mixed-effects model fit by REML
  Data: d1 
  Log-restricted-likelihood: 320824.3
  Fixed: y ~ 1 + x1 
(Intercept)          x1 
   8.302459    8.183053 

Random effects:
 Formula: ~1 + x1 | j
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr  
(Intercept) 1.699207e+00 (Intr)
x1          1.952189e+00 0.752 
Residual    1.347943e-15       

Number of Observations: 10000
Number of Groups: 100 

Answer 2

d1 <- data.marco(j, y=g00+u0j+(g10+u1j)*x1+e, x1) # Usted necesita para agregar términos de error

Mi intuición es que cuando nos simulatated un modelo mixto conjunto de datos sin los términos de error, a veces puede ser difícil para converger. Es probable que esto sea un cero residual problema.

Answer 3

Yo podría ser la falta de una respuesta más simple (es decir, errores de especificación en el modelo), pero /si/ yo estaba seguro de que el modelo mixto deberían converger, me gustaría pasar a usar openbugs, que te permitirá arrastrar mucho más mcmc de la muestra.