Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Fitting the same models in nlme and lme4

The data are from here

library(nlme)
dat0 <- read.table("aids.dat2",head=T)
dat1 <- dat0[dat0$day<=90, ]   # use only first 90-day data
dat2 <- dat1[!apply(is.na(dat1),1,any),]  # remove missing data 

# Next, let's treat the data as longitudinal (or grouped) data 
aids.dat <- groupedData(lgcopy ~ day | patid, data=dat2)

# A NLME model fit, with random effects on all 4 parameters 
start <- c(10,0.5,6,0.005)  # starting value 

aids.dat$log10copy = log10(aids.dat$lgcopy)

nlme.fit <- nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day + 1),
                 fixed = list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1),
                 random = list(patid = pdDiag(list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1))),
                 data =aids.dat, start=c(start)) 
summary(nlme.fit)

Here I fit a nonlinear mixed effects model using nlme in the nlme package. The model has 4 fixed effects and 4 random effects. I specified a diagonal structure on the variance-covariance matrix, and each patid forms a group.

library(lme4)
deriv_mod <- deriv( ~ exp(p1 - b1*t) + exp(p2 - b2*t + 1), 
                    c("p1", "b1", "p2", "b2"), function(t, p1, b1, p2, b2){})
nlmer.fit <- nlmer(deriv_mod ~ list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1) + 
                     list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1), data = aids.dat, start = c(start))

Here, I would like to fit the same model using the lme4 package. From the documentation it seems that the formula for nlmer must also have a gradient component, thus I used the deriv function first. However, I am not sure how to specify the rest of the parameters? The

deriv_mod ~ list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1) + 
                     list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1)

is to specify 4 fixed effects (in the first list object) and their corresponding 4 random effects (in the second list object). However, I am not quite sure how to specify a diagonal variance-covariance structure and make sure that the observations are grouped by patid, like I had specified in random = list(patid = pdDiag(list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1))) with nlme.

like image 909
Adrian Avatar asked Jun 26 '17 00:06

Adrian


People also ask

What is difference between LMER and Glmer?

The lmer() function is for linear mixed models and the glmer() function is for generalized mixed models.

What does lme4 do in R?

The lme4-package contains functions for estimation of multilevel or hierarchical regression models. The mlmRev-package contains, amongst many other things, the data we are going to use here. In the output below, we see that R-Project automatically loads the Matrix- and the lattice-packages as well.

What is the REML criterion?

In statistics, the restricted (or residual, or reduced) maximum likelihood (REML) approach is a particular form of maximum likelihood estimation that does not base estimates on a maximum likelihood fit of all the information, but instead uses a likelihood function calculated from a transformed set of data, so that ...

How do you know if a random effect is significant?

To determine whether a random term significantly affects the response, compare the p-value for the term in the Variance Components table to your significance level. Usually, a significance level (denoted as α or alpha) of 0.05 works well.


1 Answers

The standard way to specify fixed effects FE1 ... FE4 and independent random effects RE1 ... RE4 is shown below

mod_fit <- lme4::nlmer(Y ~ FE1 + FE2 + FE3 + FE4 + 
  (1|RE1) + (1|RE2) + (1|RE3) + (1|RE4), data= dat)

The nlme package has slightly different syntax than lme4 package.

mod_fit <- nlme::nlme(Y ~ FE1 + FE2 + FE3 + FE4 + 
      (1|RE1) + (1|RE2) + (1|RE3) + (1|RE4), 
  fixed= FE1 + FE2 + FE3 + FE4 ~ Y,
  groups= 1 ~ RE1 + RE2 + RE3 + RE4,
  data= dat)

That said, I'm not sure I completely understand the nuances of your question, so it's possible your situation means this needs to be slighly modified. If you provide comments, I'm happy to revise my answer as needed

like image 108
Alex W Avatar answered Oct 07 '22 19:10

Alex W