#Example 3.1 from Cross-over Trials in Clinical Research #8 hour PEF data of Graff-Lonevig and Browaldh #Input data n1<-7 #number of patients first sequence n2<-6 #number of patients second sequence n<-n1+n2 seqn<-factor((c(rep(1,n1),rep(2,n2),rep(1,n1), rep(2,n2))),labels=c("forsal","salfor")) #sequences patient<-factor(rep(c("1","4","6","7","0","11", "14","2","3","5","9","12","13"),2)) period<-factor(c(rep("1",n),rep("2",n))) treat<-factor((c(rep(2,n1),rep(1,n2),rep(1,n1), rep(2,n2))), labels=c("salbutamol","formoterol")) #Note: "formoterol" is coded second level of factor pef<-c(310, 310,370,410,250,380,330,370,310,380,290,260, 90,270,260,300,390,210,350,365,385,400,410,320,340,220) base<-c(290,300,250,390,250,365,190,350,350,350,280,270, 220,270,270,210,390,240,380,260,345,370,360,290,310,220) #Define contrasts options(contrasts=c(factor="contr.treatment", ordered="contr.poly")) #Demonstrate various approaches using patient as fixed effect #Analysis just fitting patient and treat fit1<-lm(pef~patient+treat) summary(fit1, corr=F) #The corr=F option suppresses the #printing out of the correlation matrix #Fitting period fit2<-lm(pef~patient+period+treat) summary(fit2,corr=F) #Fitting baselines as well fit3<-lm(pef~patient+base+period+treat) summary(fit3,corr=F) # now illustrate use of random effects model fit4<-lme(pef~seqn+period+treat, random=~1 | patient) summary(fit4) #Calculate period difference, basic estimators, #re-code seqn to seqn2 #Initialise variables pefdiff<-numeric(n) basicest<-numeric(n) seqn2<-c("forsal","salfor") #Note use of "ifelse" statement to check which sequence #patient is in # forsal sequence coded "B", salfor "A" period1 <- 1:n period2 <- (n+1):(2*n) pefdiff <- pef[period2] - pef[period1] basicest <- ifelse(seqn[period1]=="forsal", pef[period1] - pef[period2], pef[period2] - pef[period1]) seqn2 <- ifelse(seqn[period1]=="forsal","B","A") #Reset definition of contrasts options(contrasts=c(factor="contr.sum", ordered="contr.poly")) # Illustrate analysis using period differences #(equivalent to fit2 and fit4) # The regression coefficient for seqn2 will give the treatment #effect fit5<-lm(pefdiff~seqn2) summary(fit5, corr=F) #Illustrate analysis using basic estimators #(equivalent to fit2, fit4 and fit5) # The regression coefficient for intercept gives treatment #effect fit6<-lm(basicest~seqn2) summary(fit6,corr=F)