################################################## # Example on Clotting time of blood in McCullagh and Nelder (1989) # Hurn et al (1945) published data on clotting time of blood # giving clotting times in seconds for normal plasma diluted # to nine different percentage concentrations with prothrombin-free plasma # # # time = clotting time (in seconds) # conc= percentage concentration of normal plasma # lot= lot of thromboplastin (clotting agent) # # write the data > lot<-c(rep(1,9),rep(2,9)) > conc<-rep(c(5,10,15,20,30,40,60,80,100),2) > time<-c(118,58,42,35,27,25,21,19,18,69,35,26,21,18,16,13,12,12) # data matrix > data.frame(lot,conc,time) lot conc time 1 1 5 118 2 1 10 58 3 1 15 42 4 1 20 35 5 1 30 27 6 1 40 25 7 1 60 21 8 1 80 19 9 1 100 18 10 2 5 69 11 2 10 35 12 2 15 26 13 2 20 21 14 2 30 18 15 2 40 16 16 2 60 13 17 2 80 12 18 2 100 12 # some plots to choose the model # > plot(conc,time) > plot(conc,log(time)) > plot(log(conc),log(time)) > plot(log(conc/20),log(time)) > # Gamma model with log link (log-linear model) # explaining variables lot2 and log(conc/20) # # the dispersion parameter in the summary is # the estimated (coefficient of variation)**2 # = Var(res)/(E(res))**2 # > lot2<-lot-1 > fit1<-glm(time~lot2+log(conc/20),family=Gamma(link="log")) > summary(fit1) Call: glm(formula = time ~ lot2 + log(conc/20), family = Gamma(link = "log")) Deviance Residuals: Min 1Q Median 3Q Max -0.17470 -0.11596 -0.04281 0.06919 0.27749 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.69481 0.05155 71.68 < 2e-16 *** lot2 -0.47034 0.07095 -6.63 8.02e-06 *** log(conc/20) -0.58476 0.03772 -15.50 1.22e-10 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for Gamma family taken to be 0.02265044) Null deviance: 7.70867 on 17 degrees of freedom Residual deviance: 0.32110 on 15 degrees of freedom AIC: 104.28 Number of Fisher Scoring iterations: 5 > RR.estimates(fit1) RR.EST RR.LOWER RR.UPPER (Intercept) 40.2378 36.3712 44.5154 lot2 0.6248 0.5437 0.7180 log(conc/20) 0.5572 0.5175 0.6000 > # Try the interaction term # > fit2<-glm(time~lot2*log(conc/20),family=Gamma(link="log")) > summary(fit2) Call: glm(formula = time ~ lot2 * log(conc/20), family = Gamma(link = "log")) Deviance Residuals: Min 1Q Median 3Q Max -0.16943 -0.11534 -0.05440 0.07319 0.24588 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.70005 0.05416 68.317 < 2e-16 *** lot2 -0.48117 0.07659 -6.282 2.01e-05 *** log(conc/20) -0.60192 0.05462 -11.020 2.77e-08 *** lot2:log(conc/20) 0.03448 0.07724 0.446 0.662 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for Gamma family taken to be 0.02375256) Null deviance: 7.70867 on 17 degrees of freedom Residual deviance: 0.31576 on 14 degrees of freedom AIC: 105.97 Number of Fisher Scoring iterations: 5 > RR.estimates(fit2) RR.EST RR.LOWER RR.UPPER (Intercept) 40.4491 36.3754 44.9791 lot2 0.6181 0.5319 0.7182 log(conc/20) 0.5478 0.4922 0.6097 lot2:log(conc/20) 1.0351 0.8897 1.2043 > # compare the observed values to fitted values > cbind(time, fit1$fitted.values, fit2$fitted.values) time 1 118 90.509776 93.17492 2 58 60.348262 61.39093 3 42 47.609475 48.09634 4 35 40.237782 40.44915 5 27 31.744074 31.68963 6 25 26.828928 26.65106 7 21 21.165666 20.87961 8 19 17.888445 17.55981 9 18 15.700134 15.35282 10 69 56.549315 54.89945 11 35 37.704798 37.04701 12 26 29.745772 29.43289 13 21 25.140036 24.99991 14 18 19.833279 19.86178 15 16 16.762360 16.87033 16 13 13.224029 13.40304 17 12 11.176465 11.38436 18 12 9.809237 10.03041