_______________________________________________
Bendix Carstensen
Senior Statistician
Steno Diabetes Center
Niels Steensens Vej 2-4
DK-2820 Gentofte
Denmark
+45 44 43 87 38 (direct)
+45 30 75 87 38 (mobile)
bxc at steno.dk http://www.biostat.ku.dk/~bxc
www.steno.dk
###################################################
### Define some events and person-years and do the ###
traditional RR calculation
###################################################
D0 <- 15 ; D1 <- 28
Y0 <- 5532 ; Y1 <- 4783
RR <- (D1/Y1)/(D0/Y0)
erf <- exp( 1.96 * sqrt(1/D0+1/D1) )
c( RR, RR/erf, RR*erf )
exp( c( log(RR), sqrt(1/D0+1/D1) ) %*% ci.mat() )
###################################################
### Stack the data and fit a model that gives the RR
###################################################
D <- c(D0,D1) ; Y <- c(Y0,Y1); xpos <- 0:1 mm <- glm( D ~
factor(xpos), offset=log(Y), family=poisson ) exp( coef( mm ) )
###################################################
### The ci.lin function is in the Epi package
###################################################
library( Epi )
ci.lin( mm, E=T)[,5:7]
###################################################
### The same model can be fitted using the observed ### rate
as response and the person-years as weight
###################################################
mr <- glm( D/Y ~ factor(xpos), weight=Y,
family=poisson )
### Note the results are exactly the same summary(mm)$coef
summary(mr)$coef
###################################################
### This specification gives the possibility of an ###
identity link to give rate and rate difference
###################################################
ma <- glm( D/Y ~ factor(xpos), weight=Y,
family=poisson(link=identity) ) summary( ma
)$coef ci.lin( ma )[,c(1,5,6)]
###################################################
### The ci.lin function in Epi have nice facilities ### that
can be used to produce nice output
###################################################
# Scale person-years to get reasonable rates (per 1000) Y <-
Y/1000 rownames( CM ) <- c("rate 0","rate 1","RD 1 vs. 0") ma
<- glm( D/Y ~ factor(xpos), weight=Y,
family=poisson(link=identity) ) ci.lin( ma,
ctr.mat=CM )
_______________________________________________
R-sig-Epi at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-epi