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

[Ken Beath]

On 18/05/2009, at 7:35 AM, Matthew Albrecht wrote:

> 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?
>

I would interpret (1|ID:drug) as each person has a different baseline  
depending on drug, which doesn't seem physically reasonable.

I tried random effects on the shape of the polynomial (I think I've  
got the lmer commands right, but I haven't put a lot of thought into it)

lmer2<-lmer(pressure~drug*poly(timep,2)+(poly(timep,2)|ID),  
data=drugdat)
lmer3<-lmer(pressure~drug*poly(timep,2)+(poly(timep,2)|ID/drug),  
data=drugdat)

They don't improve the model fit, but that is what you would expect  
with only 18 subjects, and a residual std dev of 10.

A graph shows that there is a fair amount of error in the measurements.

xyplot(pressure~timep|ID,groups=drug,type="l",data=drugdat)

Modelling log pressure may be an alternative worth trying.

Ken

> 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")
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>