[R-sig-Epi] [epitools] Risk Difference

[BXC (Bendix Carstensen)]

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-bounces at stat.math.ethz.ch 
> [mailto:r-sig-epi-bounces at stat.math.ethz.ch] On Behalf Of BXC 
> (Bendix Carstensen)
> Sent: 26. august 2010 07:48
> To: epitools at yahoogroups.com
> Cc: r-sig-epi at stat.math.ethz.ch
> Subject: Re: [R-sig-Epi] [epitools] Risk Difference
> 
> Hi Will,
> 
> > -----Original Message-----
> > From: epitools at yahoogroups.com
> > [mailto:epitools at yahoogroups.com] On Behalf Of wleephd
> > Sent: 25. august 2010 21:23
> > To: epitools at 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)
> 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
>