CODEWITHSUNDEEP
rlme4mixed-modelsnlme
Adrian
Adrian

Reputation: 9793

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.

Upvotes: 9

Views: 2192

Answers (1)

alexwhitworth
alexwhitworth

Reputation: 4907

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

Upvotes: 3

Related Questions