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
.
The lmer() function is for linear mixed models and the glmer() function is for generalized mixed models.
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.
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 ...
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.
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With