### R code from vignette source 'sim.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: functions
###################################################
gen.data <-
function( tau = 15,
        sigma = 5,
           n1 = 50,
           n2 = 50,
           mu = 50,
         beta = -5,
        delta = 0,
        gamma = 10 )
{
mB <- mu + rep( c(0,delta), c(n1,n2) )
mF <- mB + rep( c(0,gamma), c(n1,n2) ) + beta
ai <- rnorm( n1+n2, 0, tau )
yB <- rnorm( n1+n2, mB+ai, sigma )
yF <- rnorm( n1+n2, mF+ai, sigma )
data.frame( yF, yB, bb=factor( rep(0:1,c(n1,n2)),
                               labels=c("Pl","Tr") ) )
}

wd2long <-
function( ss )
{
n1 <- table( ss$bb )[1]
n2 <- table( ss$bb )[2]
data.frame( yy = with( ss, c( yB, yF ) ),
            ii = rep( 1:(n1+n2), 2 ),
            tt = factor( rep( 1:2, each=(n1+n2) ),
                         labels=c("Bl","FU") ),
            bb = factor( rep ( rep(0:1,c(n1,n2)), 2 ),
                         labels=c("Pl","Tr") ) )
}


###################################################
### code chunk number 2: sim.rnw:44-56
###################################################
set.seed( 724368 )
n1 <- 50
n2 <- 50
ss <- gen.data( n1=n1, n2=n2, delta=0 )
ll <- wd2long( ss )

ff <- 4
SS <- gen.data( n1=ff*n1, n2=ff*n2, delta=0 )
SS <- subset( SS, yB > quantile( yB, 1-1/ff ) )
LL <- wd2long( SS )
cbind( ss=with( ss, table( bb ) ),
       SS=with( SS, table( bb ) ) )


###################################################
### code chunk number 3: sim-ex
###################################################
par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
with( ss, plot( yB, yF,
                pch=16, col=c("limegreen","red")[bb],
                xlim=c(0,100), ylim=c(0,100) ) )
with( SS, plot( yB, yF,
                pch=16, col=c("limegreen","red")[bb],
                xlim=c(0,100), ylim=c(0,100) ) )


###################################################
### code chunk number 4: fit-all
###################################################
library( lme4 )
library( Epi )
round( ci.lin( lm( yF ~ yB + bb, data = ss ) ), 4 )
round( ci.lin( lmer( yy ~ tt + I((bb=="Tr")*(tt=="FU")) + (1|ii), data=ll ) ), 4 )
round( ci.lin( lmer( yy ~ tt + bb + I((bb=="Tr")*(tt=="FU")) + (1|ii), data=ll ) ), 4 )

round( ci.lin( lm( yF ~ yB + bb, data = SS ) ), 4 )
round( ci.lin( lmer( yy ~ tt + I((bb=="Tr")*(tt=="FU")) + (1|ii), data=LL ) ), 4 )
round( ci.lin( lmer( yy ~ tt + bb + I((bb=="Tr")*(tt=="FU")) + (1|ii), data=LL ) ), 4 )