> seizure_read.table("seizure.dat") > seizure[1:10,] age base trt y1 y2 y3 y4 1 31 11 0 5 3 3 3 2 30 11 0 3 5 3 3 3 25 6 0 2 4 0 5 4 36 8 0 4 4 1 4 5 22 66 0 7 18 9 21 6 29 27 0 5 2 8 7 7 31 12 0 6 4 0 2 8 42 52 0 40 20 23 12 9 37 23 0 5 6 6 5 10 28 10 0 14 13 6 0 > names(seizure) [1] "age" "base" "trt" "y1" "y2" "y3" "y4" > model.matrix(~ age + base + trt, data=seizure)[1:10,] (Intercept) age base trt 1 1 31 11 0 2 1 30 11 0 3 1 25 6 0 4 1 36 8 0 5 1 22 66 0 6 1 29 27 0 7 1 31 12 0 8 1 42 52 0 9 1 37 23 0 10 1 28 10 0 > lmfit_lm(y4 ~ age + base + trt, data = seizure) > names(lmfit) [1] "coefficients" "residuals" "fitted.values" "effects" [5] "R" "rank" "assign" "df.residual" [9] "contrasts" "terms" "call" > summary(lmfit) Call: lm(formula = y4 ~ age + base + trt, data = seizure) Residuals: Min 1Q Median 3Q Max -15.91 -2.119 0.2008 2.028 20.23 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -4.9734 3.6234 -1.3726 0.1755 age 0.1194 0.1118 1.0678 0.2903 base 0.3078 0.0261 11.7945 0.0000 trt -1.3589 1.3743 -0.9888 0.3271 Residual standard error: 5.245 on 55 degrees of freedom Multiple R-Squared: 0.7199 F-statistic: 47.11 on 3 and 55 degrees of freedom, the p-value is 3.22e-15 Correlation of Coefficients: (Intercept) age base age -0.9369 base -0.3903 0.1884 trt -0.2871 0.0994 0.0036 > anova(lmfit) Analysis of Variance Table Response: y4 Terms added sequentially (first to last) Df Sum of Sq Mean Sq F Value Pr(F) age 1 31.910 31.910 1.1601 0.2861496 base 1 3828.834 3828.834 139.1962 0.0000000 trt 1 26.893 26.893 0.9777 0.3270968 Residuals 55 1512.871 27.507 > lmfit2_update(lmfit, ~ . - trt) > anova(lmfit2,lmfit) Analysis of Variance Table Response: y4 Terms Resid. Df RSS Test Df Sum of Sq F Value Pr(F) 1 age + base 56 1539.764 2 age + base + trt 55 1512.871 +trt 1 26.89343 0.9777029 0.3270968 > glmfit_glm(y4 ~ age + base + trt, family = "poisson", data = seizure) > names(glmfit) [1] "coefficients" "residuals" "fitted.values" [4] "effects" "R" "rank" [7] "assign" "df.residual" "weights" [10] "family" "linear.predictors" "deviance" [13] "null.deviance" "call" "iter" [16] "y" "contrasts" "terms" [19] "formula" > summary(glmfit) Call: glm(formula = y4 ~ age + base + trt, family = "poisson", data = seizure) Deviance Residuals: Min 1Q Median 3Q Max -3.163567 -1.02458 -0.1443371 0.4865469 3.899323 Coefficients: Value Std. Error t value (Intercept) 0.77557402 0.284559366 2.725526 age 0.01404377 0.008578994 1.636995 base 0.02205693 0.001088242 20.268399 trt -0.27048206 0.101852442 -2.655627 (Dispersion Parameter for Poisson family taken to be 1 ) Null Deviance: 476.2487 on 58 degrees of freedom Residual Deviance: 147.0216 on 55 degrees of freedom Number of Fisher Scoring Iterations: 4 Correlation of Coefficients: (Intercept) age base age -0.9533563 base -0.4746153 0.3092625 trt -0.3758832 0.2678822 -0.0892384 > anova(glmfit) Analysis of Deviance Table Poisson model Response: y4 Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev NULL 58 476.2487 age 1 4.4159 57 471.8328 base 1 317.7294 56 154.1033 trt 1 7.0817 55 147.0216 > glmfit2_update(glmfit, ~ . - trt) > anova(glmfit2,glmfit) Analysis of Deviance Table Response: y4 Terms Resid. Df Resid. Dev Test Df Deviance 1 age + base 56 154.1033 2 age + base + trt 55 147.0216 +trt 1 7.081701 > pvalues <- function(glmfit){ + beta <- coef(glmfit); + se <- sqrt(diag(summary(glmfit)$cov.unscaled)); + t <- beta/se; + p <- 2*(1-pnorm(abs(t))); + data.frame( + Estimate=round(beta,4), + SE=round(se,4), + t=round(t,2), + p=round(p,4)); + } > pvalues(glmfit) Estimate SE t p (Intercept) 0.7756 0.2846 2.73 0.0064 age 0.0140 0.0086 1.64 0.1016 base 0.0221 0.0011 20.27 0.0000 trt -0.2705 0.1019 -2.66 0.0079 > attach(seizure) > postscript("exploratory.plots",height=9.5,width=7,pointsize=14,horizontal=F) > par(mfrow=c(2,2)) > plot(y4 ~ base, data = seizure) > plot(log(y4+.5) ~ base, data = seizure) > hist(y4[trt==0],breaks=5*((-1):13)) > hist(y4[trt==1],breaks=5*((-1):13)) > postscript("residual.plot1",pointsize=14,horizontal=T) > par(mfrow=c(1,2)) > plot(lmfit$resid ~ lmfit$fitted) > plot(glmfit$resid ~ glmfit$fitted) > postscript("residual.plot2",pointsize=14,horizontal=T) > par(mfrow=c(1,2)) > plot(resid(lmfit) ~ predict(lmfit)) > plot(resid(glmfit) ~ predict(glmfit)) > q() # end-of-file Generated postscript file "residual.plot2". Generated postscript file "residual.plot1". Generated postscript file "exploratory.plots".