RLog to Introduce GLM Model Fitting, Model Comparison and Dispersion Extensions =============================================================================== E Slud, STAT 818R 9/25/26 Compare https://www.math.umd.edu/~slud/s705/Rlogs/RlogF09.GLM.txt Introduction of R fitting syntax: glm model-fitting function, "family=", y ~ x specification Data examples ============= install.packages("CatDataAnalysis") ## datasets for Agresti book data() ## to look at list of datasets, mostly by exercise number Data examples ============= shooting = exercise_4.13 ### 24 x 3, about basketball; same player, different games head(shooting) game made attempts 1 1 0 4 2 2 7 9 3 3 4 11 4 4 3 6 5 5 5 6 6 6 2 7 crabs = table_4.3 ### 173 x 6, about horseshoe crabs; y = indicator of "satellite" head(crabs) color spine width satell weight y 1 3 3 28.3 8 3050 1 2 4 3 22.5 0 1550 0 3 2 1 26.0 9 2300 1 4 4 3 24.8 0 2100 0 5 4 3 26.0 4 2600 1 6 3 3 23.8 0 2100 0 ## Basic fitting: shoot1 = with(shooting, made/attempts) fit1 = glm(shoot1 ~ 1, family=binomial, weight=shooting$attempts) Call: glm(formula = shoot1 ~ 1, family = binomial, weights = shooting$attempts) Coefficients: (Intercept) -0.4633 Degrees of Freedom: 23 Total (i.e. Null); 23 Residual Null Deviance: 58.31 Residual Deviance: 58.31 AIC: 96.38 ## 2nd way shoot2 = with(shooting, cbind(succ=made, fail=attempts-made)) fit2 = glm(shoot2 ~ 1, family=binomial) ### same fit, different syntax > logLik(fit1) 'log Lik.' -47.19195 (df=1) > logLik(fit2) 'log Lik.' -47.19195 (df=1) ## analysis involving dispersion, to account for erratic sequence plot(1:24,shoot1, type="b", xlab="game", ylab="field-goal rate", pch=5, main = "Ray Allen shooting across playoff games") mu = plogis(fit1$coef) ### = 0.3862069 horizontal blue line, "fitted" abline(h=mu, col="blue") ### saved as "Shooting_ex4.13.pdf" ## expected value of sum of standardized (OBS-EXP)^2 over n-p = n-1 sum( (shoot1 - mu)^2/(mu*(1-mu)/shooting$attempts) )/(24-1) ### = 2.035939 fit3 = glm(shoot2 ~ 1, family=quasibinomial) summary(fit3) ### ... same fitted coef as before, null & residual deviance etc. (Dispersion parameter for quasibinomial family taken to be 2.035939) ##==================================================================== ## Illustration of analysis of deviance for GLM using crabs data crab.fit0 = glm(y ~ color + spine + width + weight, data=crabs, family=binomial) crab.fit = glm(y ~ color + spine + width + weight, data=crabs, family=binomial(link="probit")) c(logLik(crab.fit0), logLik(crab.fit)) [1] -93.33031 -93.29981 > anova(crab.fit) Analysis of Deviance Table, Model: binomial, link: probit Response: y Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 172 225.76 color 1 12.3264 171 213.43 0.0004466 *** spine 1 1.0894 170 212.34 0.2966069 width 1 24.5338 169 187.81 7.302e-07 *** weight 1 1.2093 168 186.60 0.2714703 ## color and spine categorical (ordinal, just a few levels each) ## but width and weight continuous (many levels) crab.fitA = glm(y ~ color + spine, data=crabs, family=quasibinomial) summary(crab.fitA)$dispersion ### 1.024479 no dispersion, even though fit is horrible summary(update(crab.fitA, formula = y ~ . - spine + width)$dispersion ### = 1.00685 crab.fitB = update(crab.fitA, formula = y ~ . -spine + width*color, family=binomial) anova(crab.fitB) Analysis of Deviance Table, Model: binomial, link: logit Response: y, Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 172 225.76 color 1 12.4607 171 213.30 0.0004156 *** width 1 24.1767 170 189.12 8.789e-07 *** color:width 1 0.2221 169 188.90 0.6374498 crab.fitB = update(crab.fitA, formula = y ~ . -spine + width, family=binomial) ### So we can test for fit ... hist(crabs$width) ## nice distribution from 21 to 33.5, break at 24, 27, 29 with(crabs, table(color, cut(width, breaks=c(20.5,24,27,29,34)))) color (20.5,24] (24,27] (27,29] (29,34] 2 0 6 4 2 3 10 45 30 10 4 9 24 10 1 5 6 13 2 1 ### Make factor from cells with some pooling cell4 = with(crabs, factor(ifelse(color < 3.5, ifelse(width <= 27,"LL","LH"), ifelse(width <= 27,"HL","HH")))) num4 = table(cell4) HH HL LH LL 14 52 46 61 SSQ.fac = aggregate( (crab.fitB$y-crab.fitB$fitted)^2, by=list(cell4), sum)$x var.fac = aggregate( crab.fitB$fitted*(1-crab.fitB$fitted), by=list(cell4), sum)$x aux = cbind.data.frame(num4, SSQ.fac,var.fac, ratio=SSQ.fac/var.fac) aux cell4 Freq SSQ.fac var.fac ratio 1 HH 14 1.736660 2.417682 0.7183161 2 HL 52 10.863583 11.442920 0.9493716 3 LH 46 4.494989 4.674776 0.9615411 4 LL 61 14.820295 13.363294 1.1090300 sum(aux$ratio) ## 3.738 , compare with chisq 4df, OK !! ###--------------- install.packages("glmtoolbox") library(glmtoolbox) ### glmgee is the function we will look at Recall the Analysis of Deviance Table we produced before, likelihood-based from GLM Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 172 225.76 color 1 12.3264 171 213.43 0.0004466 *** spine 1 1.0894 170 212.34 0.2966069 width 1 24.5338 169 187.81 7.302e-07 *** weight 1 1.2093 168 186.60 0.2714703 ## We can reproduce this with a series of model fits and logLik values, as follows: LLvec = numeric(5) tmp.1 = glm(y ~ 1, data=crabs, family=binomial) LLvec[1] = logLik(tmp.1) tmp.2 = update(tmp.1, y ~ . + color); LLvec[2] = logLik(tmp.2) tmp.3 = update(tmp.2, y ~ . + spine); LLvec[3] = logLik(tmp.3) tmp.4 = update(tmp.3, y ~ . + width); LLvec[4] = logLik(tmp.4) tmp.5 = update(tmp.4, y ~ . + weight); LLvec[5] = logLik(tmp.5) -2*LLvec [1] 225.7585 213.2979 212.2519 188.0159 186.6606 ### exactly matches Res.Dev in the Analysis of deviance table ### While we are exercising our R skills this way, let's verify what is the AIC "best" ### model among all up to 2nd order interactions among the variables fit.best = step( glm(y ~ color, data=crabs, family=binomial), scope = list(lower = y ~ 1, upper = y ~ (color+spine+width+weight)^2), direction="both", k=2) Start: AIC=217.3 ### this value is: AIC(tmp.2) = 217.2979 y ~ color Df Deviance AIC + width 1 189.12 195.12 + weight 1 190.27 196.27 213.30 217.30 + spine 1 212.25 218.25 - color 1 225.76 227.76 Step: AIC=195.12y ~ color + width ### so the best final model is "color + width" ### just as we thought Df Deviance AIC 189.12 195.12 + weight 1 187.93 195.93 + spine 1 188.02 196.02 + color:width 1 188.90 196.90 - color 1 194.45 198.45 - width 1 213.30 217.30 ##============================================================================== ### NOW LET'S LOOK AT THE PURE ESTIMATING EQUATION ANALOG --- library(glmtoolbox) tmp2 = glmgee(y ~ color + spine + width + weight, data=crabs, family=binomial) > cbind(tmp$coef,tmp2$coef) ### OK, same for both functions (Intercept) -7.0078435548 -7.0078435648 color -0.5915285685 -0.5915285687 spine 0.2716837136 0.2716837134 width 0.2732524970 0.2732524974 weight 0.0007949448 0.0007949448 anova(tmp2) ### not at all the same as the glm-based anova Wald test Model 1 : y ~ 1 Model 2 : y ~ color Model 3 : y ~ color + spine Model 4 : y ~ color + spine + width Model 5 : y ~ color + spine + width + weight Chi df Pr(>Chi) 1 vs 2 11.2182 1 0.00081 *** 2 vs 3 1.0604 1 0.30312 3 vs 4 22.8225 1 1.777e-06 *** 4 vs 5 1.7281 1 0.18865 > summary(tmp.5)$coef Estimate Std. Error z value Pr(>|z|) (Intercept) -7.0078435548 3.803484545 -1.842480 0.06540499 color -0.5915285685 0.241724133 -2.447123 0.01440019 spine 0.2716837136 0.241040656 1.127128 0.25968830 width 0.2732524970 0.189330401 1.443257 0.14894798 weight 0.0007949448 0.000691723 1.149224 0.25046353 W.tsts = rbind(summary(tmp.2)$coef[2,], summary(tmp.3)$coef[3,], summary(tmp.4)$coef[4,], summary(tmp.5)$coef[5,]) W.tsts Estimate Std. Error z value Pr(>|z|) [1,] -0.7146793952 0.209473867 -3.411783 6.453945e-04 [2,] 0.2235138450 0.217958603 1.025488 3.051332e-01 [3,] 0.4572162205 0.103924124 4.399520 1.084908e-05 [4,] 0.0007949448 0.000691723 1.149224 2.504635e-01 ## Compare z^2 with chisq from glm anova table and glmgee anova : cbind(Wald=W.tsts[,3]^2, anova.chsq = anova(tmp.5)[2:5,2], anova.glmgee = anova(tmp2)[,1] ) Wald anova.chsq anova.glmgee 1 vs 2 11.640264 12.460655 11.218222 2 vs 3 1.051625 1.045974 1.060403 3 vs 4 19.355773 24.236039 22.822497 4 vs 5 1.320716 1.355241 1.728136 > ### Our likelihood theory says that the first two are ### asymptotically equivalent statistic values when ### the GLM logistic-regression model is properly specified. ### And the last column is the set of squared standardized coefficients ### using the "robustified" variance. ### Let's try to check this statement about the third column. ## I will do this for the 3 vs 4 Wald test. ## Recall that the robustified EEQ variance is B^{-1} A B^{-1} ## where B is the "observed information matrix" (negative Hessian at MLE) from likelihood th'y Binvmat = summary(tmp.4)$cov.unscaled Bmat = solve(Binvmat) rbind(sqrt(diag(Binvmat)), summary(tmp.4)$coef[,2]) (Intercept) color spine width [1,] 2.906518 0.24129 0.2404945 0.1039241 [2,] 2.906518 0.24129 0.2404945 0.1039241 ## confirms that this is how glm gets SEs ### Will discuss in class how to get A and then compute "robust" SE and Wald test for ### 3 vs 4 model in the estimating equation setting ... Amat = t(model.matrix(tmp.4)) %*% diag((crabs$y-tmp.4$fitted)^2) %*% model.matrix(tmp.4) > Amat (Intercept) color spine width (Intercept) 31.77842 110.4285 78.82632 816.9466 color 110.42846 406.3874 283.25017 2830.8212 spine 78.82632 283.2502 217.06451 2019.8577 width 816.94659 2830.8212 2019.85767 21082.7377 > Bmat (Intercept) color spine width (Intercept) 31.63706 111.5524 79.29682 816.5444 color 111.55235 413.6696 287.39507 2874.0917 spine 79.29682 287.3951 219.04106 2043.7285 width 816.54438 2874.0917 2043.72848 21168.6937 ### So these do look close. Now let's compare the Wald test chisq statistic ### produced from the likelihood assuming the correctness of the GLM model tmp.4 summary(tmp.4)$coef[4,3]^2 [1] 19.35577 ## now using the same coef estimate, get statistic standardized by `robustified' EEE sandwich variance tmp.4$coef[4]^2/(Binvmat %*% Amat %*% Binvmat)[4,4] width 22.82344 which agrees to 3 decimal places with the Wald statistic from the glmtoolbox package! ********************************************************************************************************* EXERCISE: calculate in the same way the other Wald statistics for models 1 vs 2, 2 vs 3 and 4 vs 5 using information from the tmp.1, tmp.2, tmp.3, tmp.4, tmp.5 "glm" models, to compare with the Wald statistics from the anova.glmgee output using glmtoolbox on the "crabs" data ********************************************************************************************************* ############# Now let's INTERPRET the comparison with the third column. ### Does it give evidence of fit ? of lack of fit ?