## data and analysis, after Miettinen 1976 (jh 2015 ) age.band = seq(50,75,5)+rep(c(4,9),3)/10 c= c(35,52,52,86,105,76) PT=c(77.4,68.4,61.5,47.4,38.0,23.2)*1.5 c1=c(24,35,31,46,60,39) c0=c(1,2,5,7,13,14) d1=c(22,35,38,42,51,32) d0=c(4,4,3,15,28,20) IDR.hat = (c1/d1) / (c0/d0) ; round(IDR.hat,2) V.logIDR.hat.Woolf = 1/c1+1/c0+1/d1+1/d0 ## Woolf summary ID ratio W = 1/V.logIDR.hat.Woolf log.IDR.hat = sum( W* log(IDR.hat) ) / sum(W) round( exp( log.IDR.hat + c(-1.96,0,1.96)*sqrt(1/sum(W)) ), 3) # Mantel-Haenszel n=c1+c0+d1+d0 mh.num = sum( c1*d0/n ) mh.den = sum( c0*d1/n ) round(mh.num/mh.den,3) ## crude round( ( sum(c1)*sum(d0) ) / ( sum(c0)*sum(d1) ), 2) # (prospective) logistic regression 24 rows smoke = rep( c(1,1,0,0), 6) age.stratum = rep( age.band, each = 4) Y = rep( c(1,0,1,0), 6) number = as.vector( t(cbind(c1,d1,c0,d0)) ) cbind(age.stratum,smoke,Y,number) fit=glm(Y~as.factor(age.stratum)+smoke,weights=number,family=binomial) summary(fit) cbind( round(exp(fit$coefficients),2) ) # (retrospective) logistic regression fit.retro=glm(smoke~as.factor(age.stratum)+Y,weights=number,family=binomial) cbind( round(exp(fit.retro$coefficients),2) ) ## calculate 30 year risk (paper does it differently) ID.if.smoke = (c/(c1+c0))*c1/( PT * (d1/(d0+d1)) ) ID.if.nonSm = (c/(c1+c0))*c0/( PT * (d0/(d0+d1)) ) cbind(age.band,c,PT,c1,c0,d1,d0,log.IDR.hat,W, 100*ID.if.smoke, 100*ID.if.nonSm) Thirty.year.integral.if.smoke = sum((ID.if.smoke/1000)*5) Thirty.year.integral.if.nonSm = sum((ID.if.nonSm/1000)*5) Thirty.year.risk.if.smoke = 1-exp(-Thirty.year.integral.if.smoke) Thirty.year.risk.if.nonSm = 1-exp(-Thirty.year.integral.if.nonSm) round(100*c(Thirty.year.risk.if.smoke,Thirty.year.risk.if.nonSm),2) # Risk ratio: round(Thirty.year.risk.if.smoke/Thirty.year.risk.if.nonSm,2) # this is the proper use of Risk Ratio or Relative Risk # ie a fived time horizon (30 years) ## and note that there is no followup in this study base, ## rather there is a dynamic population over 18 calendar months # jh 2015.02.12