3
$\begingroup$

I'm trying to replicate the output of mg4.math.nlme model fit via nlme() function, using the lme() function in mg5.math.lme model.

I'm wondering how I should define the fixed and random part of my lme() call to replicate the nlme() results?

ps. essentially, we want to do a multiple-group analysis using nb_wght which is a binary variable (like gender). In other words, we want means, covariances, and residual variances of the model all be specific for each value of nb_wght (i.e., the binary group).

library(nlme)

lmectl <- lmeControl(maxIter = 200,
                     msMaxIter = 200,
                     niterEM = 50,
                     msMaxEval = 400)

math <- read.csv('https://raw.githubusercontent.com/rnorouzian/v/main/mlgrp1.csv')

##========================== `nlme()` model: ==========================
mg4.math.nlme <-
  nlme(math ~ nb_wght * (
    (beta_1_N + d_1i_N) +
      (beta_2_N + d_2i_N) * (grade)
  ) +
    lb_wght * (
      (beta_1_L + d_1i_L) +
        (beta_2_L + d_2i_L) * 
        (grade)
    ),
  data = math,
  fixed = beta_1_N + beta_2_N + beta_1_L + beta_2_L ~ 1,
  random = pdBlocked(list(d_1i_N + d_2i_N ~ 1,
                          d_1i_L + d_2i_L~1)),
  group =~ id,
  start = c(35, 4, 35, 4),
  weights = varIdent(form =~ 1|factor(lb_wght)),
  control = (list(returnObject = TRUE)))

(summary(mg4.math.nlme))
#             Value  Std.Error   DF  t-value       p-value
#beta_1_N 35.481308 0.36474271  931 97.27764  0.000000e+00
#beta_2_N  4.297310 0.09017385 1287 47.65583 1.782340e-286
#beta_1_L 32.797810 1.40625509 1287 23.32280 1.232624e-100
#beta_2_L  4.878055 0.33822372 1287 14.42257  8.007625e-44

##========================== `lme()` equivalent model: ==========================
mg5.math.lme <- lme(math ~ nb_wght * grade,
                    data = math,
                    random = 
                      list(id = pdBlocked(list(
                        ~grade, ~0 + grade:nb_wght))), 
                    weights = varIdent(form =~ 1 | factor(lb_wght)),
                    control = lmectl)

(summary(mg5.math.lme))
#                   Value Std.Error   DF   t-value       p-value
#(Intercept)   32.6765937 1.3319270 1287 24.533322 2.346406e-109
#nb_wght        2.8052003 1.3815070  930  2.030536  4.258615e-02
#grade          4.9372007 0.3543493 1287 13.933146  3.332060e-41
#nb_wght:grade -0.6402842 0.3656091 1287 -1.751281  8.013572e-02
Improve this question
$\endgroup$
5
  • $\begingroup$ I don't know enough about nlme to answer, but I reformatted the code to make it easier (for me) to read, and added the AIC, which is very similar for the two models - not sure if that means something. :) $\endgroup$
    – Jeremy Miles
    Commented Feb 4, 2021 at 4:48
  • $\begingroup$ This feels like a homework question and looking at the repo with the data, it seems like a homework question. It should be tagged self-study. $\endgroup$
    – Livius
    Commented Feb 5, 2021 at 1:43
  • $\begingroup$ Even if this isn't homework, I still would like to see some explanation about how and why you converted from nlme to lme. How: what steps did you take to reformulate the problem. Why: why do you need to rewrite an nlme model as an lme model? Without some motivation, this seems like an XY Problem. If the motivation is understanding the conversion, then some of the derivation in your proposed conversion needs to be shown. $\endgroup$
    – Livius
    Commented Feb 5, 2021 at 15:15
  • $\begingroup$ Related to my previous comment on the XY Problem: what type of equivalence are you going for? Equivalent predictions? Then test the fitted values and predicted values on a hold-out set. Identical model matrices? Test those directly, although floating point error and differences in estimation between the functions may lead to some differences. One-to-one correspondence between coefficients in two different ways of expressing a model formula? Convert the first formula to math notation, do some algebra to get it into a form that lends itself to conversion to Wilkinson notation for linear models. $\endgroup$
    – Livius
    Commented Feb 5, 2021 at 15:22
  • $\begingroup$ In all these cases, it would also be helpful to know something about the structure of the data -- what are the properties of each of the data variables? Are some of them bounded? Are some categorical? Etc. $\endgroup$
    – Livius
    Commented Feb 5, 2021 at 15:24

1 Answer 1

Reset to default
3
$\begingroup$

Since this question is about re-writing a model specification, I'll focus on that aspect. Without knowing more about the data or precise inferential goal, I can't comment on whether either model is the right tool. I'm also unsure why the OP needs to rewrite a functioning model to use a different function within the same software package.

Let's start from the OP's setup.

library(nlme)
math <- read.csv('https://raw.githubusercontent.com/rnorouzian/v/main/mlgrp1.csv')

Now I'm going to take the OP's nonlinear model specification, remove some extra parentheses and use some whitespace for readability.

nonlinear <- nlme(
  math ~ nb_wght * ( (beta_1_N + d_1i_N) + (beta_2_N + d_2i_N) * grade) +
         lb_wght * ( (beta_1_L + d_1i_L) + (beta_2_L + d_2i_L) * grade),
  data = math,
  fixed = beta_1_N + beta_2_N + beta_1_L + beta_2_L ~ 1,
  random = pdBlocked(list(d_1i_N + d_2i_N ~ 1,
                          d_1i_L + d_2i_L ~ 1)),
  group = ~ id,
  start = c(35, 4, 35, 4),
  weights = varIdent(form =~ 1|factor(lb_wght)),
  control = list(returnObject = TRUE))
summary(nonlinear)

Next, I'll rename variables to more meaningful things so that I can better understand what's going on.

nonlinear <- nlme(
  math ~ nb_wght * ( (beta_nb_wght + blup_nb_wght) + (beta_nb_wght_grade + blup_nb_wght_grade) * grade) +
         lb_wght * ( (beta_lb_wght + blup_lb_wght) + (beta_lb_wght_grade + blup_lb_wght_grade) * grade),
  data = math,
  fixed = beta_nb_wght + beta_nb_wght_grade + beta_lb_wght + beta_lb_wght_grade ~ 1,
  random = pdBlocked(list(blup_nb_wght + blup_nb_wght_grade ~ 1,
                          blup_lb_wght + blup_lb_wght_grade ~ 1)),
  group = ~ id,
  start = c(35, 4, 35, 4),
  weights = varIdent(form =~ 1|factor(lb_wght)),
  control = list(returnObject = TRUE))

summary(nonlinear)

Okay, now I can see how the pieces fit together a bit better. Time to re-write it as something closer to a classical linear model by re-ordering terms.

nonlinear <- nlme(
  math ~ (beta_nb_wght + blup_nb_wght) * nb_wght + (beta_nb_wght_grade + blup_nb_wght_grade) * I(grade*nb_wght) +
         (beta_lb_wght + blup_lb_wght) * lb_wght + (beta_lb_wght_grade + blup_lb_wght_grade) * I(grade*lb_wght),
  data = math,
  # sort the coefficients in order of interaction-level:
  #  - this will make comparisons of the fixed easier later by matching sorting
  #  - note that the predictor for a two-way interaction of continuous predictors 
  #    is just the (element-wise) product of the two predictors
  fixed = beta_nb_wght + beta_lb_wght + beta_nb_wght_grade + beta_lb_wght_grade ~ 1,
  # we use the list syntax for `random` instead of depending 
  # on the separate group argument, again bringing us closer to the usual way of 
  # specifying the linear model
  random = list(id = pdBlocked(
                        list(blup_nb_wght + blup_nb_wght_grade ~ 1,
                             blup_lb_wght + blup_lb_wght_grade ~ 1))),
  start = c(35, 4, 35, 4),
  weights = varIdent(form =~ 1|factor(lb_wght)),
  control = list(returnObject = TRUE))

summary(nonlinear)

Using the linear-model formula syntax, we now have fixed effects

math ~ 0 + nb_wght + lb_wght + grade:nb_wght + grade:lb_wght # Note the lack of intercept.

and random effects:

list(~ 0 + nb_wght + grade:nb_wght, ~ 0 + lb_wght + grade:lb_wght)

The big change in both cases is to make the betas/BLUPs variables implicit and only specify the data columns.

That brings us to the linear model:

linear <- lme(
  math ~ 0 + nb_wght + lb_wght  + grade:nb_wght + grade:lb_wght,
  data = math,
  random = list(id = pdBlocked(
                        list(~ 0 + nb_wght + grade:nb_wght,
                             ~ 0 + lb_wght + grade:lb_wght))),
  weights = varIdent(form =~ 1|factor(lb_wght)),
  method = "ML", # default for nlme, but lme and lmer both default to REML
  control = list(returnObject = TRUE))

Let's compare them now.

They make make similar predictions:

> all.equal(fitted(linear), fitted(nonlinear))
[1] TRUE

The estimates of the fixed-effects are close (this comparison is why were-ordered the coefficients in the nonlinear model):

> all.equal(unname(fixed.effects(linear)), 
+           unname(fixed.effects(nonlinear)))
[1] TRUE

as are the variance parameters:

> all.equal(unname(as.numeric(VarCorr(linear)[,"StdDev"])), 
+           unname(as.numeric(VarCorr(nonlinear)[,"StdDev"])))
[1] TRUE

So these models are identical, and their log likelihoods reflect that as well:

> logLik(nonlinear)
'log Lik.' -7963.061 (df=12)
> logLik(linear)
'log Lik.' -7963.061 (df=12)

In other cases, you might find that one or the other method converges to a slightly better optimum (=better log likelihood).

If you look at the full model summaries, you'll see that they're the same model:

Non linear

> summary(nonlinear)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: math ~ (beta_nb_wght + blup_nb_wght) * nb_wght + (beta_nb_wght_grade +      blup_nb_wght_grade) * I(grade * nb_wght) + (beta_lb_wght +      blup_lb_wght) * lb_wght + (beta_lb_wght_grade + blup_lb_wght_grade) *      I(grade * lb_wght) 
  Data: math 
       AIC      BIC    logLik
  15950.12 16018.59 -7963.061

Random effects:
 Composite Structure: Blocked

 Block 1: blup_nb_wght, blup_nb_wght_grade
 Formula: list(blup_nb_wght ~ 1, blup_nb_wght_grade ~ 1)
 Level: id
 Structure: General positive-definite
                   StdDev    Corr  
blup_nb_wght       7.8921980 blp_n_
blup_nb_wght_grade 0.8795756 -0.009

 Block 2: blup_lb_wght, blup_lb_wght_grade
 Formula: list(blup_lb_wght ~ 1, blup_lb_wght_grade ~ 1)
 Level: id
 Structure: General positive-definite
                   StdDev     Corr  
blup_lb_wght       8.99415931 blp_l_
blup_lb_wght_grade 0.02853591 0.986 
Residual           5.93062152       

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | factor(lb_wght) 
 Parameter estimates:
       0        1 
1.000000 1.170217 
Fixed effects:  beta_nb_wght + beta_lb_wght + beta_nb_wght_grade + beta_lb_wght_grade ~      1 
                      Value Std.Error   DF  t-value p-value
beta_nb_wght       35.48131 0.3647427  931 97.27765       0
beta_lb_wght       32.79780 1.4062538 1287 23.32282       0
beta_nb_wght_grade  4.29731 0.0901739 1287 47.65582       0
beta_lb_wght_grade  4.87806 0.3382250 1287 14.42252       0
 Correlation: 
                   bt_nb_ bt_lb_ bt_n__
beta_lb_wght        0.000              
beta_nb_wght_grade -0.527  0.000       
beta_lb_wght_grade  0.000 -0.549  0.000

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.080667594 -0.533467953 -0.002345929  0.537432306  2.566342070 

Number of Observations: 2221
Number of Groups: 932

Linear

> summary(linear)
Linear mixed-effects model fit by maximum likelihood
  Data: math 
       AIC      BIC    logLik
  15950.12 16018.59 -7963.061

Random effects:
 Composite Structure: Blocked

 Block 1: nb_wght, nb_wght:grade
 Formula: ~0 + nb_wght + grade:nb_wght | id
 Structure: General positive-definite
              StdDev    Corr  
nb_wght       7.8921980 nb_wgh
nb_wght:grade 0.8795756 -0.009

 Block 2: lb_wght, lb_wght:grade
 Formula: ~0 + lb_wght + grade:lb_wght | id
 Structure: General positive-definite
              StdDev     Corr  
lb_wght       8.99415931 lb_wgh
lb_wght:grade 0.02853591 0.986 
Residual      5.93062152       

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | factor(lb_wght) 
 Parameter estimates:
       0        1 
1.000000 1.170217 
Fixed effects:  math ~ 0 + nb_wght + lb_wght + grade:nb_wght + grade:lb_wght 
                 Value Std.Error   DF  t-value p-value
nb_wght       35.48131 0.3647427  930 97.27765       0
lb_wght       32.79780 1.4062538  930 23.32282       0
nb_wght:grade  4.29731 0.0901739 1288 47.65582       0
lb_wght:grade  4.87806 0.3382250 1288 14.42252       0
 Correlation: 
              nb_wgh lb_wgh nb_wg:
lb_wght        0.000              
nb_wght:grade -0.527  0.000       
lb_wght:grade  0.000 -0.549  0.000

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.080667594 -0.533467953 -0.002345929  0.537432305  2.566342070 

Number of Observations: 2221
Number of Groups: 932
Improve this answer
$\endgroup$
1
  • $\begingroup$ @rnorouzian This follow-up is exactly why I asked the questions in my comments to the original question. My answer provides the algebra and that will work no matter how many predictors you have. The other bit you need is an understanding of how categorical variables are expanded into a set of numeric contrasts, then it's just repeating that algebra with those numeric contrasts. If nb_wght is just 0's and 1's, it's already equivalent to a two-level categorical variable coded with dummy/treatment contrasts. $\endgroup$
    – Livius
    Commented Feb 6, 2021 at 11:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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