#Analysis of Example 4.3 Cross-over Trials in Clinical Research
# This example cannot be analysed by SPlus or R unless the
#Venables and Ripley MASS library is used
#This analysis not illustrated in book
#This version has SPlus commands that need to be modified commented out and their R versions added
library(MASS) #Enable use of the Venables and Ripley MASS library
#Input data values
#n1 is number of values in first sequence,
#n2 is number in second sequence, n is total
#sequence is factor of sequence labels
n1<-12
n2<-12
n<-n1+n2
sequence<-factor(c(rep("forsal",n1),rep("salfor",n2)))
patient<-factor(c(3, 4, 7, 8, 9, 11, 15, 16, 19, 20, 22, 23,
1, 2, 5, 6, 10, 12, 13, 14, 17, 18, 21, 24))
EffOrd1<-c(4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 3,
4, 4, 4, 4, 4, 4, 3, 2, 2, 3)
EffOrd2<-c(4, 1, 1, 3, 4, 3, 3, 1, 3, 1, 3, 2, 4, 4,
4, 4, 4, 4, 4, 3, 4, 4, 4, 4)
#Calculate binary categories
#EffBin1<-cut(EffOrd1,breaks=c(-0.5,3.5,4.5)) # AS: in S-Plus
EffBin1<-cut(EffOrd1,breaks=c(-0.5,3.5,4.5), labels = c(0,1)) # AS: in R
#EffBin2<-cut(EffOrd2,breaks=c(-0.5,3.5,4.5)) # AS: in S-Plus
EffBin2<-cut(EffOrd2,breaks=c(-0.5,3.5,4.5), labels = c(0,1)) #AS: in R
#Calculate change scores
#ChBin<-factor(cut((EffBin1-EffBin2),breaks=c(-1.5,-0.5,0.5,1.5)),levels=c(1,2,3)) # AS: in S-Plus
ChBin<-factor(cut((as.numeric(as.character(EffBin1))-as.numeric(as.character(EffBin2))),breaks=c(-1.5,-0.5,0.5,1.5),labels=c(1,2,3))) # AS: in R
#ChOrd<-factor(cut((EffOrd1-EffOrd2),breaks=c(-3.5,-2.5,-0.5,0.5,2.5,3.5)),levels=c(1,2,3,4,5)) # AS: in S-Plus
ChOrd<-factor(cut((EffOrd1-EffOrd2),breaks=c(-3.5,-2.5,-0.5,0.5,2.5,3.5),labels=c(1,2,3,4,5))) # AS: in R
Change.frame<-data.frame(sequence,patient,EffOrd1,EffOrd2,ChBin,ChOrd)
Change.frame
#Carry out logistic regression for categorical variables
#on the change score using approach of Venabales and Ripley
#Note that estimates and SEs appear to be half what they would
#be with SAS paramaterisation
# Use proportional odds model on first change score
#fit1<-polr(ChBin~sequence) # AS: in S-Plus
fit1<-polr(ChBin~sequence, Hess = TRUE) # AS: in R (Hess = TRUE to be able to call summary())
summary(fit1)
# Use proportional odds model on second change score
## ChOrd should be ordered:
## Change.frame$ChOrdt<-ordered(Change.frame$ChOrd, levels=2:5)
#fit2<-polr(ChOrd~sequence) # AS: in S-Plus
fit2<-polr(ChOrd~sequence, Hess = TRUE) # AS: in R
summary(fit2)