[R-sig-ME] lme/lmer for drug effects on blood pressure

[Matthew Albrecht]

Dear list,
I am using lme/lmer to model blood pressure whilst under either drug or
placebo. Design- Every participant receives both drug and placebo on two
separate occasions. Blood pressure is measured in each person at 0, 75,
130, 180 and 270 min post ingestion. I have been reading Pinheiro and
Bates and the R-lists and just wanted to make sure my method is sound.
Looking at the data, it looks parabolic, so I fitted the data to a
second order polynomial:

lmer(pressure~drug*poly(timep,2)+(1|ID), data=drugdat)
or
lme(pressure~drug*poly(timep,2), random=~1|ID, data=drugdat,
na.action=na.omit)

More complicated random effects terms such as "(1|ID) + (1|ID:drug)" -
which to me means that each person has a different blood pressure
baseline, and each person's blood pressure reacts to the drug
differently(?) - make no improvements on the model (needs more
replicates?), anything I've missed or any errors in my thinking/process?

Thanks,
Matthew Albrecht
UWA, Pharmacology


Code below..............................................................

# Data generation
ID<-1:18
drug<-rep(c(1,2), c(90,90))
timep<-rep(c(0,75,130,180,270), c(18,18,18,18,18))
pressure<-c(104,128,117,115,122,122,114,107,124,88,125,97,138,126,131,133,140,111,
106,124,119,116,111,144,118,117,119,87,136,103,113,120,102,124,141,130,
92,119,117,107,107,133,108,114,114,NA,132,104,111,107,104,114,139,116,
107,126,114,96,115,142,114,120,132,103,132,102,113,122,114,123,131,115,
116,115,126,106,120,142,108,117,128,80,140,101,116,107,119,114,138,139,
104,127,128,118,140,138,120,106,126,95,118,97,107,134,123,107,136,103,
108,138,133,118,108,191,110,131,129,108,144,127,123,131,151,155,131,152,
120,143,136,123,134,164,150,130,135,NA,141,114,143,167,155,159,142,147,
135,157,141,129,141,153,136,129,149,130,NA,110,141,163,157,169,171,138,
143,153,138,129,144,160,135,130,124,114,122,114,133,140,144,153,166,132)
drugdat<-data.frame(ID, drug, timep, pressure)
drugdat$ID<-factor(drugdat$ID)
drugdat$drug<-factor(drugdat$drug)

# Quick look
with(drugdat[drugdat$pressure!="NA",], interaction.plot(timep, drug,
pressure))

# lmer/lme commands
lmer1<-lmer(pressure~drug*poly(timep,2)+(1|ID), data=drugdat)
plot(lmer1)
summary(lmer1)
anova(lmer1)

# Or the legacy version
lme1<-lme(pressure~drug*poly(timep,2), data=drugdat, random=~1|ID,
na.action=na.omit)
summary(lme1)
anova(lme1)
plot(lme1)


# More pictures if interested, I haven't figured out how the "fitted"
and the "predict" function interchange within the defined functions
below to use the lme/lmer fits yet, it is late - or early - at the moment...
lm1<-lm(pressure~drug*poly(timep,2), data=drugdat, na.action=na.omit)
plot(pressure~timep, data=drugdat, col=c("blue", "red")[drug], pch =
c(1,2)[drug])
fit<-function(x)
    predict(lm1, newdata=
        data.frame(timep=x, drug=rep(levels(drugdat$drug)[1], length(x))))
curve(fit, 0,270, add=TRUE, col="blue")
fit<-function(x)
    predict(lm1, newdata=
        data.frame(timep=x, drug=rep(levels(drugdat$drug)[2], length(x))))
curve(fit, 0,270, add=TRUE, col="red")