Enunciar el mismo modelo mixto de intercepción y pendiente aleatoria en lme que en lmer, y las ecuaciones de intercepción/pendiente aleatoria en lmer

Question

Tengo dos Q's

Q1: Tengo un modelo mixto que he formulado en lmer(), pero ahora quiero utilizar lme() porque necesito incorporar una estructura de correlación. Puedo ver que los siguientes modelos son los mismos modelos en términos de formulaciones de modelos:

fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2|Bank2),
        data = sub)
fit2 <- lme(log_age_1 ~ log_recruits + OW_P2, 
        random = ~1 + OW_P2 |Bank2,
        data=sub)

Por ejemplo, mis efectos aleatorios (también el resumen del modelo, etc.) son los mismos para lmer y lme, perfecto

> ranef(fit1)
$Bank2
  (Intercept)       OW_P2
1  1.91719166 -0.07318487
2 -3.21155864  0.26264195
3  0.05721066  0.10997179
4 -0.14943902  0.12576724
5  0.20282785 -0.01191211
6  0.46517610  0.12442915
7  1.46847545 -0.31010202
8 -0.16244538 -0.13981099
9 -0.58743867 -0.08780014
> ranef(fit2)
  (Intercept)       OW_P2
1  1.91719100 -0.07318477
2 -3.21155842  0.26264192
3  0.05720999  0.10997190
4 -0.14943971  0.12576735
5  0.20282779 -0.01191211
6  0.46517552  0.12442926
7  1.46847623 -0.31010215
8 -0.16244470 -0.13981111
9 -0.58743771 -0.08780028

Pero ahora quiero ajustar los siguientes modelos indicados en lmer() en lme();

fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2|Bank2) + (1 + log_recruits|Bank2),
            data = sub)
fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2+ log_recruits |Bank2),
            data = sub)

¿Cómo hacerlo? ¿Qué poner después de lme(...random= ? )

lme(log_age_1 ~ log_recruits + OW_P2, 
    random = ???????????????,
    data=sub)

He probado diferentes configuraciones con list() y pdDiag() etc. basadas en http://rpsychologist.com/r-guide-longitudinal-lme-lmer pero nunca encaja con la salida de mi lmer.

P2 Esta pregunta trata de obtener la ecuación correcta para cada relación lineal del paquete lmer. Consideremos el siguiente modelo con efectos aleatorios no correlacionados:

fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2|Bank2) + (1 + log_recruits|Bank2),
            data = sub)
> summary(fit1)
Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['lmerMod']
Formula: log_age_1 ~ log_recruits + OW_P2 + (1 + OW_P2 | Bank2) + (1 +      log_recruits | Bank2)
   Data: sub

REML criterion at convergence: 270.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.01579 -0.71391 -0.02338  0.54065  2.03553 

Random effects:
 Groups   Name         Variance Std.Dev. Corr 
 Bank2    (Intercept)  7.94090  2.8180        
          OW_P2        0.14820  0.3850   -1.00
 Bank2.1  (Intercept)  7.94797  2.8192        
          log_recruits 0.03295  0.1815   -1.00
 Residual              0.80904  0.8995        
Number of obs: 90, groups:  Bank2, 9

Fixed effects:
             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)    4.2703     1.7878 12.3750   2.389   0.0337 *  
log_recruits   0.6257     0.1095 14.4380   5.714 4.75e-05 ***
OW_P2         -0.4628     0.1922  8.5130  -2.408   0.0408 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) lg_rcr
log_recruts -0.652       
OW_P2       -0.713 -0.027
> coef(fit1)
$Bank2
  (Intercept) log_recruits       OW_P2
1    8.448863    0.3859443 -0.74818801
2   -2.708164    0.8123164  0.01390098
3    2.633726    0.4846902 -0.35098075
4    2.428635    0.5239703 -0.33697190
5    4.778026    0.6053660 -0.49744877
6    1.650567    0.4554392 -0.28382537
7   10.443363    0.7419798 -0.88442383
8    5.580971    0.8323623 -0.55229450
9    5.176478    0.7889412 -0.52466533
> ranef(fit1)
$Bank2
  (Intercept)       OW_P2 (Intercept) log_recruits
1   2.0892945 -0.28542162   3.7231490  -0.23972341
2  -3.4892188  0.47666737  -2.8988438   0.18664865
3  -0.8182740  0.11178563   2.1895257  -0.14097759
4  -0.9208193  0.12579449   1.5794658  -0.10169749
5   0.2538761 -0.03468239   0.3153068  -0.02030174
6  -1.3098534  0.17894102   2.6438223  -0.17022851
7   3.0865446 -0.42165744  -1.8064452   0.11631208
8   0.6553484 -0.08952811  -3.2101767   0.20669452
9   0.4531020 -0.06189894  -2.5358038   0.16327349
> fixef(fit1)
 (Intercept) log_recruits        OW_P2 
   4.2702740    0.6256677   -0.4627664

A continuación, hagamos una predicción basada en el modelo, y tracemos los valores predichos y la relación lineal del Banco2;

ggplot(sub,aes(x=OW_P2,y=log_age_1,colour=Bank2))+  geom_point() +
  geom_point(aes(y = predict(fit1)), col = "black") +
  geom_smooth(aes(y = predict(fit1), colour = Bank2), method = "lm") + facet_wrap(~Bank2)

La pregunta es cómo obtengo las estimaciones de los parámetros de alfa y beta para cada relación lineal basada en la salida del modelo anterior.

Answer 1

Para su primera pregunta, Q1, ¿intentó comparar la salida para:

fit.lmer <- lmer(log_age_1 ~ log_recruits + OW_P2 +  
                (1 + log_recruits + OW_P2 |Bank2),
                 data = sub)

y

fit.lme <- lme(log_age_1 ~ log_recruits + OW_P2, 
               random = ~ 1 + log_recruits + OW_P2|Bank2,
               data=sub)

Todo lo que tienes que hacer para fit.lme es especificar eso:

1) Las pendientes que cuantifican el efecto de log_recruits sobre log_age_1
(controlando el efecto de OW_P2) son diferentes para distintos niveles
del factor de agrupación Banco2;

2) Las pendientes que cuantifican el efecto de OW_P2 en log_age_1 (controlando el efecto de log_recruits) son diferentes para distintos niveles
del factor de agrupación Banco2.

Sólo tiene un factor de agrupación, Banco2, por lo que no es necesario complicar la especificación de la sintaxis de su modelo lme como lo haría si tuviera dos factores de agrupación cruzados (por ejemplo, Banco2 y Región).

Answer 2

Q2 (continuación) Lo anterior no parece funcionar. Aparecerá un mensaje de error

> fit2 <- lme(log_age_1 ~ log_recruits + OW_P2, 
+     random = ~ 1 + log_recruits + OW_P2|Bank2,
+     data=sub)
Error in lme.formula(log_age_1 ~ log_recruits + OW_P2, random = ~1 + log_recruits +  : 
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (10)

Los dos modelos ajustados en lmer()

fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2|Bank2) + (1 + log_recruits|Bank2),
            data = sub)
fit1 <-lmer(log_age_1 ~ log_recruits + OW_P2 +  (1 + OW_P2+ log_recruits |Bank2),
            data = sub)

corresponde a dos configuraciones diferentes de efectos aleatorios correlacionados y no correlacionados. Pero incluso la declaración de uno de ellos en lme() parece ser difícil Voy a proporcionar los dos resúmenes. Si se ven los efectos aleatorios, la estructura de correlación es muy diferente.

> summary(fit1)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of
  freedom [lmerMod]
Formula: log_age_1 ~ log_recruits + OW_P2 + (1 + OW_P2 | Bank2) + (1 +  
    log_recruits | Bank2)
   Data: sub

REML criterion at convergence: 270.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.01579 -0.71391 -0.02338  0.54065  2.03553 

Random effects:
 Groups   Name         Variance Std.Dev. Corr 
 Bank2    (Intercept)  7.94090  2.8180        
          OW_P2        0.14820  0.3850   -1.00
 Bank2.1  (Intercept)  7.94797  2.8192        
          log_recruits 0.03295  0.1815   -1.00
 Residual              0.80904  0.8995        
Number of obs: 90, groups:  Bank2, 9

Fixed effects:
             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)    4.2703     1.7878 12.3750   2.389   0.0337 *  
log_recruits   0.6257     0.1095 14.4380   5.714 4.75e-05 ***
OW_P2         -0.4628     0.1922  8.5130  -2.408   0.0408 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) lg_rcr
log_recruts -0.652       
OW_P2       -0.713 -0.027

y

> summary(fit1)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of
  freedom [lmerMod]
Formula: log_age_1 ~ log_recruits + OW_P2 + (1 + OW_P2 + log_recruits |      Bank2)
   Data: sub

REML criterion at convergence: 269.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.99533 -0.63880 -0.01517  0.49117  1.90277 

Random effects:
 Groups   Name         Variance Std.Dev. Corr       
 Bank2    (Intercept)  8.41590  2.9010              
          OW_P2        0.17082  0.4133   -0.51      
          log_recruits 0.03917  0.1979   -0.31 -0.66
 Residual              0.79974  0.8943              
Number of obs: 90, groups:  Bank2, 9

Fixed effects:
             Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)    4.1868     1.5217  7.8180   2.751 0.025555 *  
log_recruits   0.6432     0.1132 11.6770   5.680 0.000114 ***
OW_P2         -0.4739     0.1985  8.4590  -2.388 0.042399 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) lg_rcr
log_recruts -0.511       
OW_P2       -0.625 -0.312