Routine for PH 628
Lecture 12: Ordinal Logistic Regression
Gabriel Carrasco Escobar
UC San Diego; Universidad Peruana Cayetano Heredia (UPCH)Inquiries regarding this code contact gabriel.carrasco@upch.pe
rm(list=ls())
library(tidyverse); library(epiDisplay); library(tidyr); library(questionr); library(VGAM)
library(afex); library(car); library(lmtest); library(ordinal)Ordinal Logistic Regression
data <- data.frame(MAMMOGRM=c(0,0,0,0,1,1,1,1,2,2,2,2),
AGE=c(0,1,2,3,0,1,2,3,0,1,2,3),
SAMPSIZE=c(36,42,46,23,81,110,62,26,87,90,44,17))
data$MAMMOGRM <- factor(data$MAMMOGRM,
levels = c(2,1,0), # SAS interpret as a decreasing frequency
labels = c("every year", "every few year", "never"))
data$AGE <- factor(data$AGE,
levels = c(0,1,2,3),
labels = c("65-69", "70-74","75-79","80+"))
data$MAMMOGRM <- ordered(data$MAMMOGRM,c("every year", "every few year", "never"))
wtd.table(data$MAMMOGRM, data$AGE, weights = data$SAMPSIZE) # pp7 Notes## 65-69 70-74 75-79 80+
## every year 87 90 44 17
## every few year 81 110 62 26
## never 36 42 46 23# Response profile
wtd.table(data$MAMMOGRM, weights = data$SAMPSIZE)## every year every few year never
## 238 279 147# Score test for the Proportional Odds Assumption
fit1 <- vglm(MAMMOGRM ~ AGE, weights = SAMPSIZE, family=cumulative(parallel=T), data = data) # Model fit (Ordinal)
fit2 <- vglm(MAMMOGRM ~ AGE, weights = SAMPSIZE, family=cumulative(parallel=F), data = data) # Model fit (Polychotomous)
pchisq(deviance(fit1)-deviance(fit2),df=df.residual(fit1)-df.residual(fit2),lower.tail=FALSE)## [1] 0.4733083# Ho: Assumption is met
# Type III Analysis of Effects
test <- polr(MAMMOGRM ~ AGE, weights = SAMPSIZE, data = data, Hess=TRUE)
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(test,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: MAMMOGRM
## LR Chisq Df Pr(>Chisq)
## AGE 18.28 3 0.0003851 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model fit statistics
set_default_contrasts()
fit <- clm(MAMMOGRM ~ AGE, weights = SAMPSIZE, data = data)
summary(fit) ## formula: MAMMOGRM ~ AGE
## data: data
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 664 -698.62 1407.23 4(0) 6.92e-08 2.5e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AGE70-74 0.1485 0.1772 0.838 0.40193
## AGE75-79 0.6647 0.2024 3.284 0.00102 **
## AGE80+ 0.8634 0.2683 3.218 0.00129 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## every year|every few year -0.3092 0.1347 -2.295
## every few year|never 1.5721 0.1493 10.530-2*fit$logLik # -2LogL in SAS (Intercept and Covariates)## [1] 1397.235coeftest(fit)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## every year|every few year -0.30923 0.13475 -2.2949 0.022054 *
## every few year|never 1.57206 0.14930 10.5297 < 2.2e-16 ***
## AGE70-74 0.14852 0.17719 0.8382 0.402230
## AGE75-79 0.66470 0.20240 3.2840 0.001077 **
## AGE80+ 0.86337 0.26826 3.2184 0.001352 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(coef(fit)) # OR## every year|every few year every few year|never
## 0.7340132 4.8165768
## AGE70-74 AGE75-79
## 1.1601129 1.9439130
## AGE80+
## 2.3711495exp(confint(fit)) # 95% CI - OR## 2.5 % 97.5 %
## AGE70-74 0.819940 1.642808
## AGE75-79 1.308676 2.894737
## AGE80+ 1.403296 4.021965Example multivariate
# Response profile
tab1(dat$Q471)
## dat$Q471 :
## Frequency Percent Cum. percent
## every year 195 35.8 35.8
## every few year 234 43.0 78.9
## never 115 21.1 100.0
## Total 544 100.0 100.0# Score test for the Proportional Odds Assumption
fit1 <- vglm(Q471 ~ AGEGRP+INCOME+COLUMN13+EDUCAT+VISIT1, family=cumulative(parallel=T), data = dat) # Model fit (Ordinal)
fit2 <- vglm(Q471 ~ AGEGRP+INCOME+COLUMN13+EDUCAT+VISIT1, family=cumulative(parallel=F), data = dat) # Model fit (Polychotomous)
pchisq(deviance(fit1)-deviance(fit2),df=df.residual(fit1)-df.residual(fit2),lower.tail=FALSE)## [1] 0.4402203# Ho: Assumption is met
# Type III Analysis of Effects
test <- polr(Q471 ~ AGEGRP+INCOME+COLUMN13+EDUCAT+VISIT1, data = dat, Hess=TRUE)
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(test,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: Q471
## LR Chisq Df Pr(>Chisq)
## AGEGRP 11.0373 3 0.0115261 *
## INCOME 15.1052 3 0.0017289 **
## COLUMN13 2.9451 3 0.4001672
## EDUCAT 0.7495 2 0.6874615
## VISIT1 11.7716 1 0.0006014 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model fit statistics
set_default_contrasts()
fit <- clm(Q471 ~ AGEGRP+INCOME+COLUMN13+EDUCAT+VISIT1, data = dat)
summary(fit) ## formula: Q471 ~ AGEGRP + INCOME + COLUMN13 + EDUCAT + VISIT1
## data: dat
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 544 -552.34 1132.67 5(0) 4.75e-14 1.9e+02
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AGEGRP2 -0.067213 0.196567 -0.342 0.732401
## AGEGRP3 0.547973 0.231043 2.372 0.017705 *
## AGEGRP4 0.609807 0.323620 1.884 0.059520 .
## INCOME2 -0.335335 0.233816 -1.434 0.151520
## INCOME3 -0.725168 0.253966 -2.855 0.004299 **
## INCOME4 -1.065479 0.317623 -3.355 0.000795 ***
## COLUMN133 0.195077 0.252284 0.773 0.439378
## COLUMN134 -0.134585 0.250656 -0.537 0.591315
## COLUMN135 0.139298 0.279891 0.498 0.618704
## EDUCAT2 0.155499 0.294447 0.528 0.597426
## EDUCAT3 0.001672 0.281954 0.006 0.995268
## VISIT12 -0.864187 0.252969 -3.416 0.000635 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## every year|every few year -1.6298 0.4338 -3.757
## every few year|never 0.4066 0.4276 0.951-2*fit$logLik # -2LogL in SAS (Intercept and Covariates)## [1] 1104.672coeftest(fit)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## every year|every few year -1.6297868 0.4338125 -3.7569 0.0001911 ***
## every few year|never 0.4066493 0.4275676 0.9511 0.3419993
## AGEGRP2 -0.0672129 0.1965672 -0.3419 0.7325366
## AGEGRP3 0.5479734 0.2310429 2.3717 0.0180609 *
## AGEGRP4 0.6098075 0.3236202 1.8843 0.0600673 .
## INCOME2 -0.3353348 0.2338157 -1.4342 0.1521093
## INCOME3 -0.7251680 0.2539663 -2.8554 0.0044671 **
## INCOME4 -1.0654786 0.3176230 -3.3545 0.0008519 ***
## COLUMN133 0.1950769 0.2522838 0.7732 0.4397227
## COLUMN134 -0.1345851 0.2506560 -0.5369 0.5915403
## COLUMN135 0.1392985 0.2798912 0.4977 0.6189104
## EDUCAT2 0.1554991 0.2944470 0.5281 0.5976471
## EDUCAT3 0.0016724 0.2819536 0.0059 0.9952698
## VISIT12 -0.8641866 0.2529686 -3.4162 0.0006838 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exp(coef(fit)) # OR## every year|every few year every few year|never
## 0.1959714 1.5017774
## AGEGRP2 AGEGRP3
## 0.9349961 1.7297439
## AGEGRP4 INCOME2
## 1.8400771 0.7150986
## INCOME3 INCOME4
## 0.4842432 0.3445629
## COLUMN133 COLUMN134
## 1.2154045 0.8740785
## COLUMN135 EDUCAT2
## 1.1494671 1.1682409
## EDUCAT3 VISIT12
## 1.0016738 0.4213942exp(confint(fit)) # 95% CI - OR## 2.5 % 97.5 %
## AGEGRP2 0.6358489 1.3747201
## AGEGRP3 1.1008105 2.7249318
## AGEGRP4 0.9766323 3.4812653
## INCOME2 0.4516427 1.1303875
## INCOME3 0.2937389 0.7956162
## INCOME4 0.1838204 0.6395293
## COLUMN133 0.7416673 1.9959614
## COLUMN134 0.5347301 1.4299344
## COLUMN135 0.6641801 1.9918177
## EDUCAT2 0.6562065 2.0851077
## EDUCAT3 0.5764638 1.7445632
## VISIT12 0.2557609 0.6905820