Notes on polytomous response models ----------------------------------- We will make use of a pneumoconiosis data set from McCullagh and Nelder (1989). Recall the ordinal response is the severity of pneumoconiosis in coalface workers (categorized as 1 = normal, 2 = mild pneumoconiosis, 3 = severe pneumoconiosis) and the covariate is the number of years of exposure. In fact, in the following, we use log(exposure time). library(VGAM) # see hw3 p1 for how to install this data(pneumo) pneumo exposure.time normal mild severe 1 5.8 98 0 0 2 15.0 51 2 1 3 21.5 34 6 3 4 27.5 35 5 8 5 33.5 32 10 9 6 39.5 23 7 8 7 46.0 12 6 10 8 51.5 4 2 5 pneumo$let <- log(pneumo$exposure.time) ===================================================================== graphics.off() ===================================================================== # Baseline category logit model fit.bclm <- vglm(cbind(normal,mild,severe) ~ let, family=multinomial, data=pneumo) summary(fit.bclm) windows() par(mfrow=c(1,1)) matplot(pneumo$let,predict(fit.bclm,type="response"),type="l") legend(locator(),legend=c("normal","mild","severe"),lty=1:3) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.bclm,newdata=newdata,type="response") par(mfrow=c(2,1)) matplot(newdata$let, pred, type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) # Left click then center click to use locator() cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) ===================================================================== # Cumulative logit model fit.clm <- vglm(cbind(normal,mild,severe) ~ let, family=cumulative(parallel=F), data=pneumo) summary(fit.clm) windows() par(mfrow=c(1,1)) matplot(pneumo$let,predict(fit.clm,type="response"),type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.clm,newdata=newdata,type="response") par(mfrow=c(2,1)) matplot(newdata$let,pred, type="l", col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) ===================================================================== # Cumulative probit model fit.probit <- vglm(cbind(normal,mild,severe) ~ let, family=cumulative(link=probit,parallel=F), data=pneumo) summary(fit.probit) windows() par(mfrow=c(1,1)) matplot(pneumo$let,predict(fit.probit,type="response"),type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.probit,newdata=newdata,type="response") par(mfrow=c(2,1)) matplot(newdata$let, pred, type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) ===================================================================== # Proportional Odds Model fit.pom <- vglm(cbind(normal,mild,severe) ~ let, family=cumulative(parallel=T), data=pneumo) summary(fit.pom) windows() par(mfrow=c(1,1)) matplot(pneumo$let,predict(fit.pom,type="response"),type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.pom,newdata=newdata,type="response") par(mfrow=c(2,1)) matplot(newdata$let, pred, type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) ====================================================================== # Parallel vs non-parallel logits par(mfrow=c(2,1)) pred.clm <- predict(fit.clm,newdata=newdata,type="link") pred.pom <- predict(fit.pom,newdata=newdata,type="link") matplot(newdata$let,pred.clm,type="l") matplot(newdata$let,pred.pom,type="l") ===================================================================== # Comparing the fit 1-pchisq(deviance(fit.pom) - deviance(fit.clm), df=df.residual(fit.pom) - df.residual(fit.clm)) ===================================================================== # Adjacent Category Logit Model fit.aclm <- vglm(cbind(normal,mild,severe) ~ let, family=acat, data=pneumo) summary(fit.aclm) windows() par(mfrow=c(1,1)) matplot(pneumo$let,predict(fit.aclm,type="response"),type="l") # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.aclm,newdata=newdata,type="response"), par(mfrow=c(2,1)) matplot(newdata$let, pred, type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) ================================================================== windows() par(mfrow=c(2,1)) fit.aclmp <- vglm(cbind(normal,mild,severe) ~ let, family=acat(parallel=T), data=pneumo) newlet <- seq(0,8,by=.1) newdata <- cbind(pneumo[rep(1,length(newlet)),1:4],let=newlet) pred <- predict(fit.aclmp,newdata=newdata,type="response") matplot(newdata$let, pred, type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3) cumplot <- pred cumplot[,2] <- cumplot[,2] + cumplot[,3] cumplot[,1] <- cumplot[,1] + cumplot[,2] matplot(newdata$let,cumplot,type="l",col=1) # legend(locator(),legend=c("normal","mild","severe"),lty=1:3)