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

###################################################
### code chunk number 1: regr.rnw:2-6
###################################################
options( width=90,
#        prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )


###################################################
### code chunk number 2: wt-ht
###################################################
library( Epi )
data( diet )
head( diet[,c("weight","height","energy","fat","job")] )
with( diet, plot( weight ~ height, pch=16, cex=0.8,
                  xlab="Height (cm)", ylab="Weight (kg)" ) )


###################################################
### code chunk number 3: regr.rnw:175-178
###################################################
round( pnorm(1)-pnorm(-1), 3 )
round( pnorm(2)-pnorm(-2), 3 )
round( qnorm(c(0.025,0.975)), 3 )


###################################################
### code chunk number 4: regr.rnw:271-274
###################################################
library( Epi )
data( diet )
lm( weight ~ height, data=diet )


###################################################
### code chunk number 5: regr.rnw:289-291
###################################################
m1 <- lm( weight ~ height, data = diet )
summary( m1 )


###################################################
### code chunk number 6: regr.rnw:345-347
###################################################
( var(diet$weight,na.rm=TRUE) - summary(m1)$sigma^2 ) / 
  var(diet$weight,na.rm=TRUE)


###################################################
### code chunk number 7: regr.rnw:372-373
###################################################
cor( diet$weight, diet$height, use="complete.obs" )^2


###################################################
### code chunk number 8: regr.rnw:388-391
###################################################
m3 <- lm( weight ~ height + energy + fibre, data = diet )
summary( m3 )
round( ci.lin( m3 ), 3 )


###################################################
### code chunk number 9: regr.rnw:455-460
###################################################
m4 <- lm( weight ~ I(height/sd(height,na.rm=TRUE)) + 
                   I(energy/sd(energy,na.rm=TRUE)) + 
                   I(fibre /sd(fibre ,na.rm=TRUE)),
          data = diet )
summary(m4)


###################################################
### code chunk number 10: predint
###################################################
nd <- data.frame( height = c(151:190,NA,150:189),
                   fibre = rep(c(0.6,2.5),c(40,41)),
                  energy = 28 )
pr.m <- predict( m3, interval="conf", newdata=nd )
pr.p <- predict( m3, interval="pred", newdata=nd )
par( mfrow=c(1,2) )
matshade( nd$height, pr.m, plot=TRUE, ylim=c(40,100),
          xlab="Height (cm)", ylab="Weight (kg)" )
with( diet, points( height, weight, cex=0.3, pch=16 ) )
matshade( nd$height, pr.p, plot=TRUE, ylim=c(40,100),
          xlab="Height (cm)", ylab="Weight (kg)" )
with( diet, points( height, weight, cex=0.3, pch=16 ) )


###################################################
### code chunk number 11: lm-diag
###################################################
par( mfrow=c(1,3), mar=c(3,3,2,1) )
plot( m3, which=1:3, pch=16, cex=0.6, lwd=3 )


###################################################
### code chunk number 12: regr.rnw:707-709
###################################################
str( diet$job )
levels( diet$job )


###################################################
### code chunk number 13: regr.rnw:725-727
###################################################
mf <- lm( weight ~ job, data=diet )
round( ci.lin( mf ), 2 )


###################################################
### code chunk number 14: regr.rnw:759-761
###################################################
mF <- lm( weight ~ 0 + job, data=diet )
round( ci.lin( mF ), 2 )


###################################################
### code chunk number 15: regr.rnw:775-777
###################################################
mm <- lm( weight ~ job + height, data=diet )
round( ci.lin( mm ), 3 )


###################################################
### code chunk number 16: regr.rnw:793-796
###################################################
mM <- lm( weight ~ job + I(height-170), data=diet )
round( ci.lin( mM ), 3 )
anova( mM )


###################################################
### code chunk number 17: regr.rnw:833-835
###################################################
mi <- lm( weight ~ job*height, data=diet )
round( ci.lin( mi ), 3 )


###################################################
### code chunk number 18: regr.rnw:840-842
###################################################
mI <- lm( weight ~ 0 + job + job:I(height-170), data=diet )
round( ci.lin( mI ), 3 )


###################################################
### code chunk number 19: regr.rnw:852-853
###################################################
anova( mm, mi, test="Chisq" )


###################################################
### code chunk number 20: reg-int
###################################################
clr <- c("red","blue","forestgreen")
with( diet, plot( weight ~ height, 
                  pch=16, cex=0.6, col=clr[job],
                  xlab="Height (cm)", ylab="Weight (kg)" ) )
pr1 <- data.frame( height=150:190, job=levels(diet$job)[1] )
pr2 <- data.frame( height=150:190, job=levels(diet$job)[2] )
pr3 <- data.frame( height=150:190, job=levels(diet$job)[3] )
matshade( pr1$height, cbind( predict(mm,pr1,i="c"), predict(mi,pr1,i="c"),
                             predict(mm,pr2,i="c"), predict(mi,pr2,i="c"),
                             predict(mm,pr3,i="c"), predict(mi,pr3,i="c") ),
          col=rep(clr,each=2), lty=c("solid","21"), lend="butt",
          lwd=3, alpha=0.1 )
text( 155, 105-0:2*3, levels(diet$job), col=clr, font=2, adj=0 )


###################################################
### code chunk number 21: regr.rnw:936-941
###################################################
mx <- lm( weight ~ 0 + job + 
               Relevel(job,list(CB=2:3)):I(height-170), 
          data=diet )
round( ci.lin( mx )[,1:2], 3 )
anova( mm, mx, mi, test="Chisq" )


###################################################
### code chunk number 22: ggpairs
###################################################
library(ggplot2)
library(GGally)
ggpairs( diet, columns=9:13 )


###################################################
### code chunk number 23: regr.rnw:1069-1072
###################################################
mc <- lm( weight ~ height + fat + energy, data = diet )
round( cbind( ci.lin(mc)[,c(1:2,4)], 
              cor(vcov(mc)) ), 3 )


###################################################
### code chunk number 24: regr.rnw:1100-1102
###################################################
round( ci.lin( update(mc,.~.-fat   ) ), 4 )
round( ci.lin( update(mc,.~.-energy) ), 4 )


###################################################
### code chunk number 25: regr.rnw:1148-1149
###################################################
summary( update(mc,.~.+I(fat+energy)) )


###################################################
### code chunk number 26: regr.rnw:1177-1178
###################################################
summary( update(mc,.~.-fat+I(fat+energy)+fat) )


###################################################
### code chunk number 27: regr.rnw:1230-1232
###################################################
x <- c(1,1.5,2,3,9:11)
cbind(x,log2=log2(x),log10=log10(x),log1.5=log(x,base=1.5))


###################################################
### code chunk number 28: log-look
###################################################
options(digits=22)
exp(1)
options(digits=7)
curve( log, from=0.2, to=4, lwd=4, xaxs="i", xlim=c(0,4) )
abline( -1, 1 )
abline( v=1, h=0 )
abline( v=exp(1), h=1, lty="24" )
axis( side=1, at=exp(1), labels="e" )


###################################################
### code chunk number 29: regr.rnw:1318-1321
###################################################
ml <- lm( log(weight) ~ height, data=diet )
round( ci.lin( ml ), 4 )
round( (ci.exp( ml, subset="height" )-1)*100, 1 )


###################################################
### code chunk number 30: regr.rnw:1366-1367
###################################################
summary(ml)


###################################################
### code chunk number 31: regr.rnw:1404-1406
###################################################
me <- lm( weight ~ log(height,base=1.1), data=diet )
round( ci.lin( me ), 4 )


###################################################
### code chunk number 32: regr.rnw:1430-1432
###################################################
mb <- lm( log(weight) ~ log(height,base=1.1), data=diet )
round( (ci.exp( mb, subset="height" )-1)*100, 1 )


###################################################
### code chunk number 33: regr.rnw:1475-1484
###################################################
source("~/stat/R/lib.src/Epi/pkg/R/harm.R")
source("~/stat/R/lib.src/Epi/pkg/R/ci.lin.R")
class( diet$dob )
bd <- julian( diet$dob ) 
harm( bd, ord=2, per=365.25 )[1:3,]
md <- lm( height ~ harm( bd, ord=2, per=365.25 ), data=diet )
round( ci.lin(md), 3 )
Wald( md, subset=2:5 )
anova( md )


###################################################
### code chunk number 34: period
###################################################
nd <- data.frame( bd=-30:395 )
matshade( nd$bd, ci.pred(md,nd), plot=TRUE, lwd=3, las=1, 
          xaxt="n", ylim=c(170,180), 
          xlab="Season of birth", ylab="Height(cm)" )
mth.lgd <- c(31,28.25,31,30,31,30,31,31,30,31,30,31)
mth.brd <- cumsum( c(0,mth.lgd) )
mth.mid <- mth.brd[-1] - mth.lgd/2
mth.nam <- months( as.Date("2000-01-15")+0:11*30, abb=TRUE )
axis( side=1, at=mth.brd, labels=NA )
axis( side=1, at=mth.mid, labels=mth.nam, tcl=0, las=2 )
abline( v=c(0,365.25) )
abline( v=mth.brd, lty="14")

bm <- bd %% 365.25
points( bm, diet$height, pch=16, cex=0.5 )
mth <- cut( bm, mth.brd )
mmn <- tapply( diet$height, mth, mean, na.rm=TRUE )
points( mth.mid, mmn, pch=16, col="red" )


###################################################
### code chunk number 35: regr.rnw:1656-1685
###################################################
library( MASS )
N <- 10000
res <- NArray( list(sample=1:N,
                      what=c("max","min","rng","tmax","tmin","tdif") ) )
pts <- 0:365
mod.mat <- model.matrix(md,data=data.frame(bd=pts,height=175))
sim.par <- mvrnorm( n=N, coef(md), vcov(md) )

for( i in 1:N) 
   { 
sim.crv <- mod.mat %*% cbind(sim.par[i,])
res[i,"max"] <- max(sim.crv)
res[i,"min"] <- min(sim.crv)
res[i,"rng"] <- res[i,"max"] - res[i,"min"]
tmax <- pts[sim.crv==max(sim.crv)]
tmin <- pts[sim.crv==min(sim.crv)]
                       tdif <- tmax - tmin
#if( tdif > ( 365.25/2) ) tdif <- tmax -(tmin+365.25) else
#if( tdif < (-365.25/2) ) tdif <- tmax -(tmin-365.25)  
res[i,"tmax"] <- tmax
res[i,"tmin"] <- tmin
res[i,"tdif"] <- tdif
   }
ci.res <- apply( res, 2, quantile, probs=c(500,25,975)/1000 )
str( ci.res )
round( t(ci.res), 2 )
zz <- format( as.Date( t(ci.res)[4:5,], origin="2000-01-01" ), "%d %b" )
dim(zz)<-c(2,3)
zz