Incidentally a colleague of mine pointed out that you might be after a RISK difference, i.e. a difference of proportions.
Same trick: D <- c(15,28) # Numerator N <- c(22,61) # Denominator # Hand calculation: D/N diff(D/N) # Glm to give the sd easily: summary( glm( cbind(D,N-D)~factor(0:1),family=binomial(link="identity") ) ) Bendix > -----Original Message----- > From: r-sig-epi-boun...@stat.math.ethz.ch > [mailto:r-sig-epi-boun...@stat.math.ethz.ch] On Behalf Of BXC > (Bendix Carstensen) > Sent: 26. august 2010 07:48 > To: epito...@yahoogroups.com > Cc: r-sig-epi@stat.math.ethz.ch > Subject: Re: [R-sig-Epi] [epitools] Risk Difference > > Hi Will, > > > -----Original Message----- > > From: epito...@yahoogroups.com > > [mailto:epito...@yahoogroups.com] On Behalf Of wleephd > > Sent: 25. august 2010 21:23 > > To: epito...@yahoogroups.com > > Subject: [epitools] Risk Difference > > > > > > > > Hello, > > > > I am very new to R and I am trying to find a way to calculate risk > > differece and the associated CI. What is method used to > calculate CI? > > > > Thanks for you help. > > > > Will > > You can compute rate-ratios using a Poisson model with log-link. > Rate differences comes from using a Poisson model with identity link. > > Here is an illustration of how it goes; it requires small > trick in specification of the glm. The following R-code also > produces confidence intervals etc. using the function > ci.lin() from the Epi package, www.biostat.ku.dk/~bxc/Epi. > > Best regards, > Bendix Carstensen > _______________________________________________ > > 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) > b...@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@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-sig-epi > _______________________________________________ R-sig-Epi@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-epi
