### R code from vignette source 'Vickers-ex'
### Encoding: UTF-8

###################################################
### code chunk number 1: acup-read
###################################################
library( Epi )
library( foreign )
acp <- read.dta( # "c:/bendix/artikler/rep-meas/data/sportsmen.dta" )[,-4]
                 "./data/sportsmen.dta" )[,-4]
names( acp ) <- c("bl","fu","gr")
acp$gr <- factor( acp$gr, labels=c("Placebo","Acupuncture") )
str( acp )
head( acp )


###################################################
### code chunk number 2: Vickers-ex.rnw:58-60
###################################################
mc <- lm( fu ~ bl + gr, data=acp )
summary( mc )


###################################################
### code chunk number 3: Vickers-ex.rnw:70-71
###################################################
( mb <- mean( acp$bl ) )


###################################################
### code chunk number 4: Vickers-ex.rnw:76-78
###################################################
( cf <- coef( mc ) )
(cf-c(0,1,1)) %*% cbind( c(1,mb,0),  c(1,mb,1) )


###################################################
### code chunk number 5: Vickers-ex.rnw:83-85
###################################################
(cf-c(0,1,1)) %*% cbind( c(1,40,0),  c(1,40,1) )
(cf-c(0,1,1)) %*% cbind( c(1,75,0),  c(1,75,1) )


###################################################
### code chunk number 6: acup-raw
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
with( acp, plot( bl, fu, pch=16, col=c("blue","red")[gr],
                 xlim=c(20,100), ylim=c(20,100),
                 xlab="Baseline score", ylab="Follow-up score" ) )
abline( 0, 1 )
text( rep(100,2), c(25,30), levels(acp$g), font=2, col=c("blue","red"), adj=1 )


###################################################
### code chunk number 7: acup-reg
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
with( acp, plot( bl, fu, pch=16, col=c("blue","red")[gr],
                 xlim=c(20,100), ylim=c(20,100),
                 xlab="Baseline score", ylab="Follow-up score" ) )
cf <- coef( mc )
abline( 0, 1 )
abline( v=mb )
abline( v=c(40,75), lty="26" )
abline( cf[1]      , cf[2], lwd=3, col="blue" )
abline( cf[1]+cf[3], cf[2], lwd=3, col="red" )
text( rep(100,2), c(25,30), levels(acp$g), font=2, col=c("blue","red"), adj=1 )


###################################################
### code chunk number 8: acup-fu
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
with( acp, plot( bl, fu, pch=16, col=c("blue","red")[gr],
                 xlim=c(20,100), ylim=c(20,100),
                 xlab="Baseline score", ylab="Follow-up score" ) )
( fu <- with( acp, tapply( fu   , gr, mean ) ) )
diff( fu )
abline( 0, 1 )
abline( v=mb )
abline( v=c(40,75), lty="26" )
abline( h=fu[1], lwd=2, lty=2, col="blue" )
abline( h=fu[2], lwd=2, lty=2, col="red" )
text( rep(100,2), c(25,30), levels(acp$g), font=2, col=c("blue","red"), adj=1 )


###################################################
### code chunk number 9: acup-dif
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
with( acp, plot( bl, fu, pch=16, col=c("blue","red")[gr],
                 xlim=c(20,100), ylim=c(20,100),
                 xlab="Baseline score", ylab="Follow-up score" ) )
( df <- with( acp, tapply( fu-bl, gr, mean ) ) )
diff( df )
abline( 0, 1 )
abline( v=mb )
abline( v=c(40,75), lty="26" )
abline( df[1], 1, lwd=2, lty=2, col="blue" )
abline( df[2], 1, lwd=2, lty=2, col="red" )
text( rep(100,2), c(25,30), levels(acp$g), font=2, col=c("blue","red"), adj=1 )


###################################################
### code chunk number 10: acup
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
with( acp, plot( bl, fu, pch=16, col=c("blue","red")[gr],
                 xlim=c(20,100), ylim=c(20,100),
                 xlab="Baseline score", ylab="Follow-up score" ) )
df <- with( acp, tapply( fu-bl, gr, mean ) )
fu <- with( acp, tapply( fu   , gr, mean ) )
abline( 0, 1 )
abline( v=mb )
abline( v=c(40,75), lty="26" )
abline( cf[1]      , cf[2], lwd=3, col="blue" )
abline( cf[1]+cf[3], cf[2], lwd=3, col="red" )
abline( h=fu[1], lwd=2, lty=2, col="blue" )
abline( h=fu[2], lwd=2, lty=2, col="red" )
abline( df[1], 1, lwd=2, lty=2, col="blue" )
abline( df[2], 1, lwd=2, lty=2, col="red" )
text( rep(100,2), c(25,30), levels(acp$g), font=2, col=c("blue","red"), adj=1 )


###################################################
### code chunk number 11: Vickers-ex.rnw:186-202
###################################################
# Follow-up
fu <- with( acp, tapply( fu   , gr, mean ) )
c( fu, diff( fu ) )
mf <- lm( fu ~ gr, data=acp )

# differences
df <- with( acp, tapply( fu-bl, gr, mean ) )
c( df, diff( df ) )
md <- lm( fu-bl ~ gr, data=acp )

# compare with regression
cmp.cf <- rbind( ci.lin( mc, subset="Acu" ),
                 ci.lin( mf, subset="Acu" ),
                 ci.lin( md, subset="Acu" ) )
rownames( cmp.cf ) <- c("Cond.","FU onl","Diff.")
round( cmp.cf, 4 )


###################################################
### code chunk number 12: Vickers-ex.rnw:209-212
###################################################
lg <- reshape( acp, varying=1:2, v.names="score", direction="long" )
head( lg )
str( lg )


###################################################
### code chunk number 13: Vickers-ex.rnw:217-220
###################################################
library( lme4 )
mr <- lmer( score ~  gr + gr:factor(time) + (1|id), data=lg )
round( summary( mr )$coef, 2 )


###################################################
### code chunk number 14: Vickers-ex.rnw:230-232
###################################################
mR <- lmer( score ~  gr*factor(time) + (1|id), data=lg )
round( summary( mR )$coef, 2 )


###################################################
### code chunk number 15: Vickers-ex.rnw:236-241
###################################################
summary( mR )
VarCorr( mR )
( tausq <- as.numeric( VarCorr( mR )$id ) )
( sigsq <- attr( VarCorr( mR ), "sc" )^2 )
( rho <- tausq/(tausq+sigsq) )


###################################################
### code chunk number 16: Vickers-ex.rnw:244-247
###################################################
round( summary( mR )$coef, 2 )
round( cf <- fixef( mR ), 2 )
( 1- rho ) * cf[2] + cf[4]


###################################################
### code chunk number 17: Vickers-ex.rnw:253-261
###################################################
lg <- transform( lg, G2 = Relevel( interaction( gr, time ),
                                   list( B=1:2 ) ) )
with( lg, ftable( gr, time, G2 ) )
ms <- lmer( score ~ G2 + (1|id), data=lg )
round( summary( ms )$coef, 3 )
# CM <- rbind( diag(3), c(0,-1,1) )
# rownames( CM ) <- c( names( fixef(ms) ), "Acp-eff" )
# round( ci.lin( ms, ctr.mat=CM ), 2 )


###################################################
### code chunk number 18: Vickers-ex.rnw:270-272
###################################################
round( ci.lin( lm( fu    ~ gr, data=acp ) )[2,], 2 )
round( ci.lin( lm( fu-bl ~ gr, data=acp ) )[2,], 2 )