#------------------------------------------------------------------- # 8. GAMMA FAMILY AND CENSORED SURVIVAL DATA #------------------------------------------------------------------- #------------------------------------------------------------------- # INSULATE #------------------------------------------------------------------- > m <- matrix(scan("/glim/datasets/insulate"), ncol = 2, byrow = T) > kv <- m[, 1] > time <- m[, 2] > par(mfrow = c(2, 2)) > plot(kv, log(time), main = "Figure 8.1 Insulate") > summary(glm(time ~ 1, family = Gamma(link = log))) Coefficients: Value Std. Error t value (Intercept) 4.590641 0.396575 11.57572 (Dispersion Parameter for Gamma family taken to be 11.95265 ) Residual Deviance: 371.6379 on 75 degrees of freedom Number of Fisher Scoring Iterations: 0 > glm <- glm(time ~ kv, family = Gamma(link = log)) > summary(glm) Coefficients: Value Std. Error t value (Intercept) 21.4118567 1.61239375 13.27955 kv -0.5553653 0.04845234 -11.46210 (Dispersion Parameter for Gamma family taken to be 1.734156 ) Residual Deviance: 131.6428 on 74 degrees of freedom Number of Fisher Scoring Iterations: 6 > chisq <- (371.6379 - 131.6428)/1.734156 > abline(21.4118567, -0.5553653) > c(chisq, 1 - pchisq(chisq, 1)) [1] 138.393 0.000 > kv2 <- kv^2 > summary(glm(time ~ kv + kv2, family = Gamma(link = log))) Coefficients: Value Std. Error t value (Intercept) 16.548968415 14.05844068 1.1771553 kv -0.250501247 0.86724728 -0.2888464 kv2 -0.004730961 0.01329415 -0.3558679 (Dispersion Parameter for Gamma family taken to be 1.715972 ) Residual Deviance: 131.4256 on 73 degrees of freedom Number of Fisher Scoring Iterations: 6 > chisq <- (131.6428 - 131.4256)/1.715972 > c(chisq, 1 - pchisq(chisq, 1)) [1] 0.1265755 0.7220095 > u <- 1 - exp( - time/fitted(glm)) > hist(u, main = "Figure 8.2 Uniform?") > u <- sort(u) > uhat <- ((1:76) - 0.5)/76 > plot(uhat, u, xaxs = "s", yaxs = "s", main = "Figure 8.3 P-P plot") > abline(0, 1) > k <- sqrt(76) * (max(abs(u - uhat)) + 0.5/76) > k [1] 1.265117 #------------------------------------------------------------------- # USING THE INVERSE LIKELIHOOD #------------------------------------------------------------------- > ystar <- 1/time > summary(glm(ystar ~ 1, family = poisson, weight = time)) Coefficients: Value Std. Error t value (Intercept) -4.590641 0.4732592 -9.700058 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 371.6379 on 75 degrees of freedom Number of Fisher Scoring Iterations: 0 > summary(glm(ystar ~ kv, family = poisson, weight = time)) Coefficients: Value Std. Error t value (Intercept) -21.4108609 1.26254262 -16.95853 kv 0.5553353 0.03795337 14.63204 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 131.6428 on 74 degrees of freedom Number of Fisher Scoring Iterations: 5 #------------------------------------------------------------------- # CENSORED EXPONENTIAL DATA #------------------------------------------------------------------- > censor <- (time < 10) * 1 > censor.time <- censor * time + (1 - censor) * 10 > censor.time [1] 5.79 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 [12] 10.00 10.00 10.00 10.00 10.00 10.00 10.00 7.74 0.40 10.00 9.88 [23] 10.00 10.00 2.75 0.79 10.00 3.91 0.27 0.69 10.00 10.00 10.00 [34] 10.00 0.96 4.15 0.19 0.78 8.01 10.00 7.35 6.50 8.27 10.00 [45] 10.00 3.16 4.85 2.78 4.67 1.31 10.00 10.00 10.00 1.97 0.59 [56] 2.58 1.69 2.71 10.00 0.35 0.99 3.99 3.67 2.07 0.96 5.35 [67] 2.90 10.00 0.47 0.73 1.40 0.74 0.39 1.13 0.09 2.38 > plot(kv, log(time), main = "Figure 8.4 Censored data", pch = ".") > points(kv, log(censor.time)) > abline(log(10), 0, lty = 2) > abline(21.4118567, -0.5553653) > ystar <- censor/censor.time > summary(glm(ystar ~ 1, family = poisson, weight = censor.time)) Coefficients: Value Std. Error t value (Intercept) -2.353256 0.1546169 -15.21991 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 156.3604 on 75 degrees of freedom Number of Fisher Scoring Iterations: 0 > summary(glm(ystar ~ kv, family = poisson, weight = censor.time)) Coefficients: Value Std. Error t value (Intercept) -17.5702038 2.433342 -7.220607 kv 0.4556724 0.069908 6.518172 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 103.4691 on 74 degrees of freedom Number of Fisher Scoring Iterations: 5 > abline(17.5702038, -0.4556724, lty = 3) #------------------------------------------------------------------- # LEUKEMIA #------------------------------------------------------------------- > m <- matrix(scan("/glim/datasets/leukemx"), ncol = 3, byrow = T) > sample <- factor(m[, 1]) > weeks <- m[, 2] > censor <- m[, 3] > rate <- censor/weeks > summary(glm(rate ~ 1, family = poisson, weight = weeks)) Coefficients: Value Std. Error t value (Intercept) -2.892222 0.2675086 -10.8117 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 54.50253 on 41 degrees of freedom Number of Fisher Scoring Iterations: 0 > summary(glm(rate ~ sample, family = poisson, weight = weeks)) Coefficients: Value Std. Error t value (Intercept) -2.9227905 0.199101 -14.679941 sample 0.7633062 0.199101 3.833764 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 38.01732 on 40 degrees of freedom Number of Fisher Scoring Iterations: 5 > summary(glm(rate ~ 1, family = binomial(link = cloglog), weight = weeks )) Coefficients: Value Std. Error t value (Intercept) -2.863833 0.2675449 -10.70412 (Dispersion Parameter for Binomial family taken to be 1 ) Residual Deviance: 61.50737 on 41 degrees of freedom Number of Fisher Scoring Iterations: 0 > summary(glm(rate ~ sample, family = binomial(link = cloglog), weight = weeks)) Coefficients: Value Std. Error t value (Intercept) -2.8861191 0.199141 -14.49284 sample 0.7873098 0.199141 3.95353 (Dispersion Parameter for Binomial family taken to be 1 ) Residual Deviance: 43.96778 on 40 degrees of freedom Number of Fisher Scoring Iterations: 6 #------------------------------------------------------------------- # MYELOMA #------------------------------------------------------------------- > m <- matrix(scan("/glim/datasets/myeloma"), ncol = 7, byrow = T) > censor <- m[, 1] > time <- m[, 2] > urea <- m[, 3] > hemo <- m[, 4] > age <- m[, 5] > sex <- factor(m[, 6]) > calcium <- m[, 7] > rate <- censor/time > options(contrasts = c("contr.treatment", "contr.poly")) > summary(glm(rate ~ urea + hemo + age + sex + calcium, family = binomial( link = cloglog), weight = time)) Coefficients: Value Std. Error t value (Intercept) -4.56494060 1.62962137 -2.8012277 urea 1.70587065 0.56513790 3.0185034 hemo -0.13031944 0.06022672 -2.1638144 age -0.02137618 0.01552488 -1.3768984 sex -0.06087744 0.31016824 -0.1962723 calcium 0.14297425 0.09674710 1.4778143 (Dispersion Parameter for Binomial family taken to be 1 ) Residual Deviance: 72.95315 on 59 degrees of freedom Number of Fisher Scoring Iterations: 6 > summary(glm(rate ~ urea, family = binomial(link = cloglog), weight = time)) Coefficients: Value Std. Error t value (Intercept) -5.763102 0.7811063 -7.378128 urea 1.711872 0.5509226 3.107283 (Dispersion Parameter for Binomial family taken to be 1 ) Residual Deviance: 79.80026 on 63 degrees of freedom Number of Fisher Scoring Iterations: 6 #------------------------------------------------------------------- # AML #------------------------------------------------------------------- > m <- matrix(scan("/glim/datasets/aml1"), ncol = 13, byrow = T) > age <- m[, 1] > surv <- m[, 2] > censor <- m[, 3] > rig <- m[, 4] > rigg <- m[, 5] > rigm <- m[, 6] > riga <- m[, 7] > aig <- m[, 8] > aigg <- m[, 9] > aigm <- m[, 10] > aiga <- m[, 11] > wbc <- m[, 12] > lim <- m[, 13] > rate <- censor/surv > summary(glm(rate ~ 1, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -6.625329 0.1608219 -41.19668 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 62.46444 on 40 degrees of freedom Number of Fisher Scoring Iterations: 0 > summary(glm(rate ~ age + rig + rigg + rigm + riga + aig + aigg + aigm + aiga + wbc + lim, family = poisson, weight = surv)) Coefficients: (2 not defined because of singularities) Value Std. Error t value (Intercept) -8.9071373672 1.856845887 -4.79691795 age -0.0006724838 0.016807169 -0.04001172 rig -1.2054364863 1.143383837 -1.05427106 rigg 1.4164425707 1.280751062 1.10594682 rigm 1.4229622432 1.418068390 1.00345107 riga NA NA NA aig 0.0064617610 0.008175686 0.79036315 aigg -0.0091044329 0.009667850 -0.94172257 aigm -0.0075942033 0.009680756 -0.78446390 aiga NA NA NA wbc 0.3238993719 0.245423519 1.31975686 lim 0.0385391939 0.017732306 2.17338875 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 47.36696 on 31 degrees of freedom Number of Fisher Scoring Iterations: 9 > summary(glm(rate ~ age + rig + aig + wbc + lim, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -9.6652095113 1.8392826961 -5.25487981 age 0.0004803002 0.0153718552 0.03124543 rig 0.0558569601 0.1374272956 0.40644735 aig -0.0011210720 0.0007808975 -1.43561995 wbc 0.4202789316 0.2317480283 1.81351675 lim 0.0457535700 0.0180093401 2.54054672 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 49.59376 on 35 degrees of freedom Number of Fisher Scoring Iterations: 9 > summary(glm(rate ~ age + rig + wbc + lim, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -7.646859046 1.10605122 -6.9136573 age 0.005247831 0.01522182 0.3447571 rig -0.113457571 0.07229416 -1.5693878 wbc 0.096914819 0.05261260 1.8420457 lim 0.025347686 0.01073623 2.3609474 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 51.56908 on 36 degrees of freedom Number of Fisher Scoring Iterations: 9 > summary(glm(rate ~ age + aig + wbc + lim, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -9.0518415954 1.0191878239 -8.8814264 age 0.0015543385 0.0151038121 0.1029103 aig -0.0008504435 0.0004104394 -2.0720317 wbc 0.3443562727 0.1379062098 2.4970324 lim 0.0402556657 0.0118025295 3.4107660 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 49.75637 on 36 degrees of freedom Number of Fisher Scoring Iterations: 9 > summary(glm(rate ~ age + wbc + lim, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -8.724161765 0.95118386 -9.1718985 age 0.008195491 0.01577934 0.5193812 wbc 0.082410752 0.05322539 1.5483355 lim 0.029811295 0.01043439 2.8570220 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 54.52516 on 37 degrees of freedom Number of Fisher Scoring Iterations: 8 > summary(glm(rate ~ wbc + lim, family = poisson, weight = surv)) Coefficients: Value Std. Error t value (Intercept) -8.41148798 0.70998795 -11.847367 wbc 0.08858353 0.05085871 1.741757 lim 0.02951225 0.01027507 2.872220 (Dispersion Parameter for Poisson family taken to be 1 ) Residual Deviance: 54.7898 on 38 degrees of freedom Number of Fisher Scoring Iterations: 8