R-code/Ordinal Regression.R at master · apculley/R-code · GitHub

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
air<-airquality

# taken from https://stats.idre.ucla.edu/r/dae/ordinal-logistic-regression/

library(ordinal)

str(air)
#add a categorical variable to the data
air$Temp_cat<-air$Temp<85
air$Temp_cat<-factor(air$Temp_cat,c('TRUE','FALSE'))
levels(air$Temp_cat)<-c('Cold','Hot')
# set month as an ordinal factor
air$Month<-factor(air$Month,ordered = T)

### Method 1, using ordinal package ###
colnames(air)
fit<-clm(Month~Temp_cat+Ozone+Solar.R+Wind,data=air, link='logit')

convergence(fit)

# nominal test = for ordered independent vars, if sig, odds not proportional bw levels
nominal_test(fit)
# scale test = for unordered independent vars, if sig, odds not proportional bw levels (this model fails, unsurprisingly)
scale_test(fit)

levels(air$Month)
sf <- function(y) {
  c('Y>=5' = qlogis(mean(y >= 1)),
    'Y>=6' = qlogis(mean(y >= 2)),
    'Y>=7' = qlogis(mean(y >= 3)),
    'Y>=8' = qlogis(mean(y >= 4)),
    'Y>=9' = qlogis(mean(y >= 5)))
}
require(Hmisc)
# if proportional, diff bw levels should be equal across levels
(s <- with(air, summary(as.numeric(Month) ~ Temp_cat+Ozone+Solar.R+Wind, fun=sf)))
# distances bw points should be same for all levels of each factor (clearly not the case for this model)
plot(s, which=1:5, pch=1:5, xlab='logit', main=' ',xlim=c(-5,5))

# check predictions from model
predict(fit, type='class')

### Method 2, using
library(MASS)
m <- polr(Month~Temp_cat+Ozone+Solar.R+Wind,data=air, Hess=TRUE)

# summary doesn't output p-values, so need to calculate them separately
summary(m)
(ctable <- coef(summary(m)))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
(ci <- confint(m))
# convert to odds ratios
exp(coef(m))
exp(cbind(OR = coef(m), ci))

levels(air$Month)
sf <- function(y) {
  c('Y>=5' = qlogis(mean(y >= 1)),
    'Y>=6' = qlogis(mean(y >= 2)),
    'Y>=7' = qlogis(mean(y >= 3)),
    'Y>=8' = qlogis(mean(y >= 4)),
    'Y>=9' = qlogis(mean(y >= 5)))
}
require(Hmisc)
# if proportional, diff bw levels should be equal across levels
(s <- with(air, summary(as.numeric(Month) ~ Temp_cat+Ozone+Solar.R+Wind, fun=sf)))
# distances bw points should be same for all levels of each factor (clearly not the case for this model)
plot(s, which=1:5, pch=1:5, xlab='logit', main=' ',xlim=c(-5,5))

# check predictions from model
predict(m, type='class')