-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathKielOrdLogReg.R
141 lines (109 loc) · 5.08 KB
/
KielOrdLogReg.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
wants <- c("MASS", "ordinal", "rms", "VGAM")
has <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])
## ------------------------------------------------------------------------
set.seed(123)
N <- 100
X1 <- rnorm(N, 175, 7)
X2 <- rnorm(N, 30, 8)
Ycont <- 0.5*X1 - 0.3*X2 + 10 + rnorm(N, 0, 6)
Yord <- cut(Ycont, breaks=quantile(Ycont), include.lowest=TRUE,
labels=c("--", "-", "+", "++"), ordered=TRUE)
dfOrd <- data.frame(X1, X2, Yord)
## ------------------------------------------------------------------------
library(VGAM)
(vglmFit <- vglm(Yord ~ X1 + X2, family=propodds, data=dfOrd))
## ----results='hide'------------------------------------------------------
vglm(Yord ~ X1 + X2, family=cumulative(parallel=TRUE, reverse=TRUE), data=dfOrd)
# not shown
## ----results='hide'------------------------------------------------------
vglm(Yord ~ X1 + X2, family=acat(parallel=TRUE), data=dfOrd)
# not shown
## ----results='hide'------------------------------------------------------
vglm(Yord ~ X1 + X2, family=sratio(parallel=TRUE), data=dfOrd)
# not shown
## ------------------------------------------------------------------------
library(rms)
(ormFit <- orm(Yord ~ X1 + X2, data=dfOrd))
## ----results='hide'------------------------------------------------------
library(MASS)
(polrFit <- polr(Yord ~ X1 + X2, method="logistic", data=dfOrd))
# not shown
## ------------------------------------------------------------------------
exp(MASS:::confint.polr(polrFit))
## ----results='hide'------------------------------------------------------
library(ordinal)
(clmFit <- clm(Yord ~ X1 + X2, link="logit", data=dfOrd))
# not shown
## ------------------------------------------------------------------------
PhatCateg <- VGAM::predict(vglmFit, type="response")
head(PhatCateg)
## ----results='hide'------------------------------------------------------
predict(ormFit, type="fitted.ind")
predict(clmFit, subset(dfOrd, select=c("X1", "X2"), type="prob"))$fit
predict(polrFit, type="probs")
# not shown
## ------------------------------------------------------------------------
categHat <- levels(dfOrd$Yord)[max.col(PhatCateg)]
head(categHat)
## ----results='hide'------------------------------------------------------
predict(clmFit, type="class")
predict(polrFit, type="class")
# not shown
## ------------------------------------------------------------------------
Nnew <- 3
dfNew <- data.frame(X1=rnorm(Nnew, 175, 7),
X2=rnorm(Nnew, 30, 8))
## ------------------------------------------------------------------------
VGAM::predict(vglmFit, dfNew, type="response")
## ----results='hide'------------------------------------------------------
predict(ormFit, dfNew, type="fitted.ind")
predict(polrFit, dfNew, type="probs")
predict(clmFit, subset(dfNew, select=c("X1", "X2"), type="prob"))$fit
# not shown
## ------------------------------------------------------------------------
facHat <- factor(categHat, levels=levels(dfOrd$Yord))
cTab <- xtabs(~ Yord + facHat, data=dfOrd)
addmargins(cTab)
## ------------------------------------------------------------------------
(CCR <- sum(diag(cTab)) / sum(cTab))
## ------------------------------------------------------------------------
VGAM::deviance(vglmFit)
VGAM::logLik(vglmFit)
VGAM::AIC(vglmFit)
## ------------------------------------------------------------------------
vglm0 <- vglm(Yord ~ 1, family=propodds, data=dfOrd)
LLf <- VGAM::logLik(vglmFit)
LL0 <- VGAM::logLik(vglm0)
## ------------------------------------------------------------------------
as.vector(1 - (LLf / LL0))
## ------------------------------------------------------------------------
as.vector(1 - exp((2/N) * (LL0 - LLf)))
## ------------------------------------------------------------------------
as.vector((1 - exp((2/N) * (LL0 - LLf))) / (1 - exp(LL0)^(2/N)))
## ------------------------------------------------------------------------
sumOrd <- summary(vglmFit)
(coefOrd <- coef(sumOrd))
## ------------------------------------------------------------------------
zCrit <- qnorm(c(0.05/2, 1 - 0.05/2))
(ciCoef <- t(apply(coefOrd, 1, function(x) x["Estimate"] - zCrit*x["Std. Error"] )))
## ----results='hide'------------------------------------------------------
summary(polrFit)
summary(clmFit)
# not shown
## ------------------------------------------------------------------------
vglmR <- vglm(Yord ~ X1, family=propodds, data=dfOrd)
VGAM::lrtest(vglmFit, vglmR)
## ------------------------------------------------------------------------
VGAM::lrtest(vglmFit, vglm0)
## ------------------------------------------------------------------------
vglmP <- vglm(Yord ~ X1 + X2, family=cumulative(parallel=TRUE, reverse=TRUE),
data=dfOrd)
# vglmNP <- vglm(Yord ~ X1 + X2, family=cumulative(parallel=FALSE, reverse=TRUE),
# data=dfOrd)
# VGAM::lrtest(vglmP, vglmNP)
## ------------------------------------------------------------------------
clmP <- clm(Yord ~ X1 + X2, link="logit", data=dfOrd)
## model with non-proportional odds for X2:
clmNP <- clm(Yord ~ X1, nominal=~X2, data=dfOrd)
anova(clmP, clmNP)