2
$\begingroup$

I have two Q's

Q1: I have a mixed model that I stated in lmer(), but now I want to use lme() because I need to incorporate a correlation structure. I can see that the following models are the same models in terms of model formulations:

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)

For example, my random effects (also model summary etc.) are the same for both lmer and lme, perfect

> 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

But now I want to fit following models stated in lmer() in 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)

How to do that? what to put in after lme(...random= ?)

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

I have tried different setups with list() and pdDiag() etc. based on http://rpsychologist.com/r-guide-longitudinal-lme-lmer , but it never fits with my lmer output.

Q2 This Q is about getting the right equation for each linear relationship from the lmer package. Let's consider following model with uncorrelated random effects:

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

Then let us predict based on the model, and plot the predicted values and linear relationship of Bank2;

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)

enter image description here

Question is how do I get the parameter estimates for alpha and beta for each linear relationship based on above model output?

$\endgroup$
1
$\begingroup$

For your first question, Q1, did you try comparing the output for:

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

and

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

All you need to do for fit.lme is to specify that:

1) The slopes quantifying the effect of log_recruits on log_age_1
(controlling for the effect of OW_P2) are different for different levels
of the grouping factor Bank2;

2) The slopes quantifying the effect of OW_P2 on log_age_1 (controlling for the effect of log_recruits) are different for different levels
of the grouping factor Bank2.

You only have one grouping factor, Bank2, so there is no need to complicate the syntax specification for your lme model the way you would if you had two crossed grouping factors (e.g., Bank2 and Region).

$\endgroup$
  • $\begingroup$ It do not seem to work: $\endgroup$ – fishy_stuff Oct 1 '18 at 8:47
  • $\begingroup$ See answer below $\endgroup$ – fishy_stuff Oct 1 '18 at 8:55
0
$\begingroup$

Q2 (continued) The above do not seem to work. Error-message will occur

> 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)

The two models fitted in 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)

corresponds to two different set-ups of correlated and uncorrelated random effects. But even stating one of them in lme() seems to be difficult I'll provide the two summaries. If you see the random effects then correlation structure are very different.

> 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

and

> 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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.