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

###################################################
### code chunk number 1: prev.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: prev.rnw:22-25
###################################################
library( Epi )
data( pr )
str( pr )


###################################################
### code chunk number 3: prev.rnw:42-43
###################################################
subset( pr, sex=="M" & A==60 )


###################################################
### code chunk number 4: prev.rnw:67-69
###################################################
set.seed( 1952 )
rbinom( 1, 100, 0.1055323 )


###################################################
### code chunk number 5: prev.rnw:165-168
###################################################
p <- 16/100
se <- sqrt(p*(1-p)/100)
round( c( p, p-2*se, p+2*se ), 4 )


###################################################
### code chunk number 6: prev.rnw:175-178
###################################################
P <- 3/25
SE <- sqrt(P*(1-P)/25)
round( c( P, P-2*SE, P+2*SE ), 4 )


###################################################
### code chunk number 7: prev.rnw:220-224
###################################################
erf <- exp(1.96/sqrt(16*84/100))
pL <- p/erf / ( 1 - p + p/erf )
pU <- p*erf / ( 1 - p + p*erf )
round( c( p, pL, pU ), 4 )


###################################################
### code chunk number 8: prev.rnw:281-284
###################################################
pnorm(q=2,mean=0,sd=1) - pnorm(q=-2,mean=0,sd=1)
c( qnorm(p=0.025), qnorm(p=0.975) )
c( qnorm(p=0.05) , qnorm(p=0.95) )


###################################################
### code chunk number 9: prev.rnw:320-327
###################################################
x <- 16
n <- 100
m0 <- glm( cbind(x,n-x) ~ 1, family=binomial(link=logit) )
round( ci.lin(m0), 4 )
1/sqrt(100*0.16*0.84)
odds <- ci.exp( m0 )
round( odds / (1+odds), 4 )


###################################################
### code chunk number 10: prev.rnw:350-352
###################################################
ml <- glm( cbind(x,n-x) ~ 1, family=binomial(link=log) )
ci.exp(ml)


###################################################
### code chunk number 11: prev.rnw:430-434
###################################################
p <- 16/100
( se.p <- sqrt(p*(1-p)/n) )
tt.p <- ( p - 0.1 ) / se.p
2*pnorm( abs(tt.p), lower=FALSE )


###################################################
### code chunk number 12: prev.rnw:454-459
###################################################
ml <- glm( cbind(x,n-x) ~ 1, family=binomial(link=log) )
ci.lin(ml)
tt.m <- (-1.832581-log(0.10)) / 0.2291288
2* pnorm( abs(tt.m), lower.tail=FALSE )
Wald( ml, H0=log(0.1) )


###################################################
### code chunk number 13: prev.rnw:523-526
###################################################
prt <- stat.table( index = list(A,sex), 
                contents = ratio(X,N,100), data = pr )[1,,]
round( prt[70:75,], 1 )


###################################################
### code chunk number 14: raw
###################################################
matplot( 0:99+0.5, prt, type="l", lwd=2, col=c(4,2), lty=1 )
abline( v=c(76,82))


###################################################
### code chunk number 15: prev.rnw:579-586
###################################################
pr$A <- pr$A+0.5
ma <- glm( cbind(X,N-X) ~ Ns( A, knots=c(10,30,50,60,70,80,90) ), 
           family = binomial, 
             data = subset(pr,sex=="M") )
fa <- glm( cbind(X,N-X) ~ Ns( A, knots=c(10,30,50,60,70,80,90) ), 
           family = binomial, 
             data = subset(pr,sex=="F") )


###################################################
### code chunk number 16: step-smooth
###################################################
matplot( 0:99, prt, yaxs="i", ylim=c(0,20), 
         xlab="Age", ylab="Prevalence of diabetes 2010 (%)",
         type="s", lwd=1, col=c("blue","red"), lty=1 )
nd <- data.frame( A=seq(0,100,0.5) )
pr.m <- ci.pred( ma, nd ) * 100
pr.f <- ci.pred( fa, nd ) * 100
matshade( nd$A, cbind(pr.m,pr.f), col=c("blue","red"), lwd=1, alpha=0.2 )


###################################################
### code chunk number 17: prev.rnw:755-764
###################################################
m1 <- glm( cbind(X,N-X) ~ Ns( A, knots=c(10,   50,      70,      90) ), 
           family = binomial, 
           data = subset(pr,sex=="M") )
m2 <- glm( cbind(X,N-X) ~ Ns( A, knots=c(10,30,50,60,   70,   80,90) ), 
           family = binomial, 
           data = subset(pr,sex=="M") )
m3 <- glm( cbind(X,N-X) ~ Ns( A, knots=c(10,30,50,60,65,70,75,80,90) ), 
           family = binomial, data = subset(pr,sex=="M") )
anova( m1, m2, m3, test="Chisq" )


###################################################
### code chunk number 18: prev.rnw:797-799
###################################################
c( m1$deviance, m1$null.deviance )       
c( nrow(m1$data), m1$df.residual, m1$df.null )


###################################################
### code chunk number 19: prev.rnw:820-822
###################################################
summary(m1)
c( m1$aic, m2$aic, m3$aic )


###################################################
### code chunk number 20: prev.rnw:842-844
###################################################
library( stats4 )
c( BIC(m1), BIC(m2), BIC(m3) )