I would like do create a mixed linear model for an unbalanced dataset (different number of events per subject and a few missing values for some time points). I am using R version 3.2.1 (2015-06-18)
, package: nlme_3.1-120
.
Here comes simulated data:
library(nlme)
set.seed(1)
subject <- factor(rep(c(1, 1, 2, 3, 4, 4, 4, 5, 6, 7, 7, 8, 9, 9, 10,
11, 11, 11, 12, 13), 10))
event <- factor(rep(1:20, 10))
timepoint <- rep(1:10, each = 20)
measure <- rnorm(length(timepoint)) + timepoint*0.3
timepoint <- factor(timepoint)
measure[sample(1:length(measure), rpois(5,4))] <- NA
data <- data.frame(subject=subject, event=event, timepoint=timepoint,
measure=measure)
str(data)
The model should predict the variable “measure” over different time points as fixed effect and for subjects and events as random effects.
base <- lme(measure ~ 1, data=data, random= ~ 1|subject,
na.action=na.exclude, method="ML")
intercept <- lme(measure ~ timepoint, data=data, random= ~ 1|subject,
na.action=na.exclude, method="ML")
nested <- lme(measure ~ timepoint, data=data, random= ~ 1|subject/event,
na.action=na.exclude, method="ML")
anova(base, intercept, nested)
I would like to fit random intercept and slope, because intercept and slope can vary among subjects and events. However when I add the random slope effect, the model does not converge. It does not through any error message, but it runs to infinity. What can I do create a model with random slope that converges?
cave model runs endless
slope <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject,
na.action=na.exclude, method="ML")
I tried also this
cave model runs endless
slope2 <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject,
na.action=na.exclude, method="ML", control=list(opt="optim"))
cave some models may run endless
slope3 <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject/event,
na.action=na.exclude, method="ML", control = list(opt="optim"))
covariance <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject,
correlation=corAR1(),na.action = na.exclude, method="ML")
covariance2 <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject,
correlation=corAR1(0), na.action=na.exclude, method="ML",
control=list(opt="optim"))
covariance3 <- lme(measure ~ timepoint, data=data, random= ~ timepoint|subject,
correlation=corAR1(0), na.action=na.exclude, method="ML",
control=list(maxlter=1000))
Best Answer
@AdamO has done a good job identifying the specific error in your code. Let me address the question more generally. Here is how I simulate a linear mixed effects model:
Mixed effects models assume each unit has random effects drawn from a multivariate normal distribution. (When a model is estimated, it is the variances and covariances of that multivariate normal that are being estimated for the random effects.) I start by specifying this distribution and generating (pseudo-)random values to serve as the random effects. It is often convenient to specify the variances as $1$, so that the covariance is the correlation between slopes and intercepts (which is easier for me to conceptualize).
Next, I would generate my $X$ variables. I can't really follow the logic of your example, so I will use
time
as my only regressor.Having generated your random effects and your regressors, you can specify the data generating process. Since you want some randomly missed timepoints, there is a level of additional complexity here. (Note that these data are missing completely at random; for more on simulating missing data, see: How to simulate the different types of missing data.)
At this point, you can fit your model. I typically use the
lme4
package.