cancer.dta file availabe at http://www.sph.emory.edu/~dkleinb/logreg2.htm.
The grade variable is an ordinal variable (natural ordering), whereas the subtype is a multinomial variable (no natural ordering).
Both have three-levels, so a regular logistic regression model cannot fit this as an outcome variable.
library(foreign)
cancer <- read.dta("http://www.sph.emory.edu/~dkleinb/datasets/cancer.dta")
cancer <- within(cancer, {
grade <- factor(grade, 0:2, c("well differentiated","moderately differentiated","poorly differentiated"),
ordered = TRUE)
race <- factor(race, 0:1, c("white","black"))
estrogen <- factor(estrogen, 0:1, c("never used","ever used"))
subtype <- factor(subtype, 0:2, c("Adenocarcinoma","Adenosquamous","Other"))
age <- factor(age, 0:1, c("50-64","65-79"))
smoking <- factor(smoking, 0:1, c("non-smoker","current smoker"))
})
Comparisons are carried out at each possible cutoff (indicated by VS below) along the ordinal scale (two for 3-level ordinal outcome variable). It is assumed that in each of the comparison below, the odds ratio assessing the effect of an exposure variable will be the same.
Comparison 1
“well differentiated” (ref) VS “moderately differentiated” and “poorly differentiated”
Comparison 2
“well differentiated” and “moderately differentiated” (ref) VS “poorly differentiated”
This function is consistent with SAS PROC LOGISTIC with “DESCENDING option” if used with “reverse = TRUE”.
Ordinal logistic regression without proportional odds assumptions (two coefficients for each predictor)
## Load VGAM
library(VGAM)
## Without proportional odds assumptions (two coefficients for each predictor)
vglm.race.estrogen.nonparallel <- vglm(grade ~ race + estrogen, cancer,
family = cumulative(link = "logit", parallel = FALSE, reverse = TRUE))
summary(vglm.race.estrogen.nonparallel)
Call:
vglm(formula = grade ~ race + estrogen, family = cumulative(link = "logit",
parallel = FALSE, reverse = TRUE), data = cancer)
Pearson Residuals:
Min 1Q Median 3Q Max
logit(P[Y>=2]) -1.583 -0.833 0.319 0.904 1.38
logit(P[Y>=3]) -0.829 -0.636 -0.277 -0.232 2.60
Coefficients:
Estimate Std. Error z value
(Intercept):1 0.585 0.231 2.53
(Intercept):2 -1.391 0.278 -5.00
raceblack:1 0.376 0.307 1.22
raceblack:2 0.516 0.351 1.47
estrogenever used:1 -0.874 0.270 -3.24
estrogenever used:2 -0.552 0.345 -1.60
Number of linear predictors: 2
Names of linear predictors: logit(P[Y>=2]), logit(P[Y>=3])
Dispersion Parameter for cumulative family: 1
Residual deviance: 574.32 on 566 degrees of freedom
Log-likelihood: -287.16 on 566 degrees of freedom
Number of iterations: 4
## Odds ratios
exp(coef(vglm.race.estrogen.nonparallel))
(Intercept):1 (Intercept):2 raceblack:1 raceblack:2 estrogenever used:1 estrogenever used:2
1.79412 0.24876 1.45608 1.67530 0.41731 0.57595
## Check parallelism
is.parallel(vglm.race.estrogen.nonparallel)
(Intercept) race estrogen
FALSE FALSE FALSE
Ordinal logistic regression with proportional odds assumptions (same coefficients for each predictor)
## With proportional odds assumptions (same coefficients for each predictor)
vglm.race.estrogen.parallel <- vglm(grade ~ race + estrogen, cancer,
family = cumulative(link = "logit", parallel = TRUE, reverse = TRUE))
summary(vglm.race.estrogen.parallel)
Call:
vglm(formula = grade ~ race + estrogen, family = cumulative(link = "logit",
parallel = TRUE, reverse = TRUE), data = cancer)
Pearson Residuals:
Min 1Q Median 3Q Max
logit(P[Y>=2]) -1.563 -0.847 0.323 0.964 1.34
logit(P[Y>=3]) -0.847 -0.585 -0.292 -0.223 2.75
Coefficients:
Estimate Std. Error z value
(Intercept):1 0.511 0.215 2.38
(Intercept):2 -1.274 0.229 -5.57
raceblack 0.427 0.272 1.57
estrogenever used -0.776 0.249 -3.11
Number of linear predictors: 2
Names of linear predictors: logit(P[Y>=2]), logit(P[Y>=3])
Dispersion Parameter for cumulative family: 1
Residual deviance: 575.21 on 568 degrees of freedom
Log-likelihood: -287.61 on 568 degrees of freedom
Number of iterations: 4
## Check parallelism
is.parallel(vglm.race.estrogen.parallel)
(Intercept) race estrogen
FALSE TRUE TRUE
Check proportional odds assumption
## Difference in -2 log likelihood
(diff.neg.2logLik <- -2*(logLik(vglm.race.estrogen.parallel) - logLik(vglm.race.estrogen.nonparallel)))
[1] 0.89729
## Difference in DF
(diff.df <- df.residual(vglm.race.estrogen.parallel) - df.residual(vglm.race.estrogen.nonparallel))
[1] 2
## Likelihood ratio test
pchisq(q = diff.neg.2logLik, df = diff.df, lower.tail = FALSE)
[1] 0.63849
## Likelihood ratio test was actually implemented in VGAM
lrtest(vglm.race.estrogen.nonparallel,
vglm.race.estrogen.parallel)
Likelihood ratio test
Model 1: grade ~ race + estrogen
Model 2: grade ~ race + estrogen
#Df LogLik Df Chisq Pr(>Chisq)
1 566 -287
2 568 -288 2 0.9 0.64
The lines marked “:1” are the ones for the first comparison, in this case levels 2 and 3 vs level 1. The lines marked “:2” are the ones for the second comparison, in this case level 3 (Other) vs levels 2 and 1.
The “parallel = FALSE” option to the cumulative() function allows different slopes for comparisons, thus, comparing this model to a proportional odds model (different intercepts, but same slopes for comparisons) tests for the proportional odds assumption.
Here the proportional odds model is adequate. The Black race compared to the White race is at 1.53 times higher odds of having a higher grade cancer at either one of the two thresholds adjusting for estrongen use status. The estrogen ever users compared to the never users are at 0.46 times lower odds of having a higher grade cancer at either one of the two thresholds adjusting for race.
One of the categories of the outcome variable is designated as the reference category, and each of the other levels is compared with this reference level. In this example the outcome level Adenocarcinoma is chosen as the reference. So the comparison will be between Adenosquamous vs Adenocarcinoma and between Other vs Adenocarcinoma.
Comparison 1
“Adenocarcinoma” (ref) VS “Adenosquamous”
Comparison 2
“Adenocarcinoma” (ref) VS “Other”
The VGAM::vglm() function can fit multinomial logistic regression as a special case of a vector generalized linear models.
One line of code for each model specification, but actually two models are constructed simultaneously. The one that compares the reference outcome level Adenocarcinoma to Adenosquamous, and the other that compares the refence outcome level Adenocarcinoma to Others.
## Load VGAM
library(VGAM)
## Using the first level as the reference level for all comparison
vglm.age.est.smk <- vglm(subtype ~ age + estrogen + smoking, data = cancer,
family = multinomial(parallel = FALSE, refLevel = 1))
summary(vglm.age.est.smk)
Call:
vglm(formula = subtype ~ age + estrogen + smoking, family = multinomial(parallel = FALSE,
refLevel = 1), data = cancer)
Pearson Residuals:
Min 1Q Median 3Q Max
log(mu[,2]/mu[,1]) -0.981 -0.472 -0.398 -0.180 3.97
log(mu[,3]/mu[,1]) -0.632 -0.592 -0.514 -0.218 5.14
Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.882 0.402 -4.676
(Intercept):2 -1.203 0.319 -3.772
age65-79:1 0.987 0.412 2.397
age65-79:2 0.282 0.328 0.861
estrogenever used:1 -0.644 0.344 -1.874
estrogenever used:2 -0.107 0.307 -0.349
smokingcurrent smoker:1 0.889 0.525 1.693
smokingcurrent smoker:2 -1.791 1.046 -1.713
Number of linear predictors: 2
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
Dispersion Parameter for multinomial family: 1
Residual deviance: 494.4 on 564 degrees of freedom
Log-likelihood: -247.2 on 564 degrees of freedom
Number of iterations: 5
The lines marked “:1” are the ones for the first comparison, in this case level 2 (Adenosquamous) vs level 1 (Adenocarcinoma). The lines marked “:2” are the ones for the second comparison, in this case level 3 (Other) vs level 1 (Adenocarcinoma).
Interpretation Age is statistially significant for the Adenosquamous vs Adenocarcinoma comparison (comparison 1), but not for the Other vs Adenoarcinoma comparison (comparison 2) after adjustment for estrogen use and smoking status. People aged 65-79 (compared to people aged 50-64) were at 2.68 times more odds of having Adenosquamous (compared to Adenocarcinoma).
Comparison 1: Adenosquamous vs Adenocarcinoma model
The estimated log odds of Adenosquamous in comparison to Adenocarcinoma among those who are aged 65-79 is as follows. The odds ratio of having Adenosquamous in comparison to Adenocarcinoma among those who are aged 65-79 compared to those are aged 50-64 is exp(0.98706) = 2.68.
ln\( [\frac {\hat{P}(D = AdenoSq|X_1)} {\hat{P}(D = AdenoCa|X_1)}] = -1.88218 + (0.98706) \)
Comparison 2: Other vs Adenocarcinoma model
The estimated log odds of Other in comparison to Adenocarcinoma among those who are aged 65-79 is as follows. The odds ratio of having Other in comparison to Adenocarcinoma among those who are aged 65-79 compared to those are aged 50-64 is exp(0.28229) = 1.33.
ln\( [\frac {\hat{P}(D = Other|X_1)} {\hat{P}(D = AdenoCa|X_1)}] = -1.20322 + (0.28229) \)
Constraining parameters
Some parameters can be assumed to be the same across the different comparisons. In this case the coefficient of the estrogen variable is assumed to be the same across the comparison, thus, only one estimate for this variable, whereas the others have two estimates, one for each comparison.
vglm.age.est_para.smk <- vglm(subtype ~ age + estrogen + smoking, data = cancer,
family = multinomial(parallel = ~ -1 + estrogen, refLevel = 1))
summary(vglm.age.est_para.smk)
Call:
vglm(formula = subtype ~ age + estrogen + smoking, family = multinomial(parallel = ~-1 +
estrogen, refLevel = 1), data = cancer)
Pearson Residuals:
Min 1Q Median 3Q Max
log(mu[,2]/mu[,1]) -0.908 -0.509 -0.377 -0.212 3.60
log(mu[,3]/mu[,1]) -0.663 -0.575 -0.487 -0.209 5.49
Coefficients:
Estimate Std. Error z value
(Intercept):1 -2.013 0.395 -5.097
(Intercept):2 -1.081 0.300 -3.598
age65-79:1 0.989 0.410 2.414
age65-79:2 0.273 0.329 0.832
estrogenever used -0.340 0.250 -1.360
smokingcurrent smoker:1 0.864 0.521 1.657
smokingcurrent smoker:2 -1.787 1.047 -1.708
Number of linear predictors: 2
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
Dispersion Parameter for multinomial family: 1
Residual deviance: 496.14 on 565 degrees of freedom
Log-likelihood: -248.07 on 565 degrees of freedom
Number of iterations: 5
When compared by the likelihood ratio test, it appears that the simpler model with single coefficient for the estrogen variable is adequate (the fit is not significantly different).
lrtest(vglm.age.est.smk, vglm.age.est_para.smk)
Likelihood ratio test
Model 1: subtype ~ age + estrogen + smoking
Model 2: subtype ~ age + estrogen + smoking
#Df LogLik Df Chisq Pr(>Chisq)
1 564 -247
2 565 -248 1 1.73 0.19