guaguncher
guaguncher

Reputation: 77

How can I produce an estimate for a changing variable using lme4/nlme?

I am fairly new to R and mixed models analysis.

I wanted to produce a single estimate for the change in the variable ln_ahr based on time for every individual. I believe this can be thought of as a slope of change in time. Here is how my data is structured (long format):

       v001   ln_ahr    time
13404 28337 0.28438718    0
13405 28337         NA    3
13406 28337         NA    6
13407 28337 1.05015991    9
13408 28337         NA   12
13409 28337 1.33345188   15
13410 28337         NA   19
13413 28355 1.14904314    0
13414 28355         NA    3
13415 28355 1.06546008    6
13416 28355         NA    9
13417 28355 1.17865500   12
13418 28355 2.84949593   15
13423 29983 0.07015499    0
13424 29983 0.21056477    3
13426 29983 0.36125306    9
13427 29983 0.66139848   12
13428 29983 0.16962391   16

where v001 is the subject identifier.

I tried to calculate the slope using the nlme package in R as:

   slope <- lme(ln_ahr~time,random=~1+time|v001,
                data=restructured,na.action="na.omit")

and I tried obtained the ranef(slope) and coef(slope) values. I read that the coef(slope) values "computes the sum of the fixed and random effects coefficients for each explanatory variable of each grouping factor" thus I believed that printing out the coefficients for time (leaving out the intercept values) would give me an estimate of each individual's change in ln_ahr over time and I can use that as my "slope" or estimate for change in ln_ahr.

Time is calculated as years where time 0 indicates the first year of ln_ahr measurement; everyone is measured every three years.

I am wondering if this a correct approach at all or if I did it correctly; if not what are your suggestions?

Upvotes: 0

Views: 594

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226087

The basic answer is "yes"; the numbers returned by lme4::coef() are the estimated subject-specific parameters.

Variants of this example are around on the internet, but:

Fit one of the basic random-slopes examples from lme4:

library(lme4)
fm1 <- lmer(Reaction~Days+(Days|Subject),sleepstudy)

Extract estimated intercepts and slopes for each group (Subject):

d1 <- coef(fm1)$Subject
d1$Subject <- rownames(d1)

For comparison, fit a separate model for each group and extract subject-specific slopes and intercepts:

fm2 <- lmList(Reaction~Days|Subject,sleepstudy)
d2 <- coef(fm2)
d2$Subject <- rownames(d2)

Plot (random-effects estimates as solid lines; fixed-effects dashed)

library(ggplot2); theme_set(theme_bw())
gg0 <- ggplot(sleepstudy,aes(Days,Reaction,
                      colour=Subject))+
  geom_point()+
  geom_abline(data=d1,
        mapping=aes(intercept=`(Intercept)`,
                    slope=Days,colour=Subject))+
  geom_abline(data=d2,linetype=2,
              mapping=aes(intercept=`(Intercept)`,
                          slope=Days,colour=Subject))  

enter image description here

Doug Bates doesn't like these "spaghetti plots" with all groups in a single panel, he prefers facets:

gg0+facet_wrap(~Subject)+
 theme(panel.spacing=grid::unit(0,"lines"))          

enter image description here

(Ideally we would also order the subjects in some non-arbitrary way, e.g. by slope)

Upvotes: 2

Related Questions