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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
library(ordinal)
#################################
## 1 categorical variable in nominal:
fm <- clm(rating ~ temp, nominal=~contact, data=wine)
fm$Theta
fm$alpha.mat
## Threshold effects:
fm <- clm(rating ~ temp, nominal=~contact, data=wine,
threshold="symmetric")
fm$Theta
fm$alpha.mat
fm <- clm(rating ~ temp, nominal=~contact, data=wine,
threshold="equidistant")
fm$Theta
fm$alpha.mat
## Singular fit is still ok (with a warning, though)
fm <- clm(rating ~ contact, nominal=~temp, data=wine)
fm$alpha.mat
fm$Theta
#################################
## 1 continuous variable:
set.seed(123)
x <- rnorm(nrow(wine), sd=1)
fm <- clm(rating ~ temp, nominal=~ x, data=wine)
fm$alpha.mat
fm$Theta
fm <- clm(rating ~ temp, nominal=~ poly(x, 2), data=wine)
fm$alpha.mat
fm$Theta
#################################
## 1 categorical + 1 continuous variable:
set.seed(123)
x <- rnorm(nrow(wine), sd=1)
fm <- clm(rating ~ temp, nominal=~contact + x, data=wine)
fm$alpha.mat
fm$Theta
fm <- clm(rating ~ temp, nominal=~contact + x, data=wine,
threshold="symmetric")
fm$alpha.mat
fm$Theta
#################################
### NOTE: To get the by-threshold nominal effects of continuous terms
## use:
with(fm, t(apply(alpha.mat, 1, function(th) tJac %*% th)))
#################################
## Interactions:
fm <- clm(rating ~ temp, nominal=~contact:x, data=wine)
fm$alpha.mat
fm$Theta
fm <- clm(rating ~ temp, nominal=~contact+x+contact:x, data=wine)
fm$alpha.mat
fm$Theta
fm <- clm(rating ~ temp, nominal=~contact*x, data=wine)
fm$alpha.mat
fm$Theta
## polynomial terms:
fm <- clm(rating ~ temp, nominal=~contact + poly(x, 2), data=wine)
fm$alpha.mat
fm$Theta
## logical variables: (treated like numeric variables)
wine$Con <- as.character(wine$contact) == "yes"
fm <- clm(rating ~ temp, nominal=~Con, data=wine)
fm$Theta
fm$alpha.mat
wine$Con.num <- 1 * wine$Con
fm <- clm(rating ~ temp, nominal=~Con.num, data=wine)
fm$Theta
fm$alpha.mat
#################################
## Two continuous variables:
set.seed(321)
y <- rnorm(nrow(wine), sd=1)
fm1 <- clm(rating ~ temp, nominal=~y + x, data=wine)
fm1$alpha.mat
fm1$Theta
## summary(fm1)
#################################
## 1 categorical + 2 continuous variables:
fm1 <- clm(rating ~ temp, nominal=~y + contact + x, data=wine)
fm1$alpha.mat
fm1$Theta
fm1 <- clm(rating ~ temp, nominal=~contact + x + contact:x + y,
data=wine)
summary(fm1)
fm1$Theta
fm1$alpha.mat
fm1 <- clm(rating ~ temp, nominal=~contact*x + y, data=wine)
fm1$Theta
fm1$alpha.mat
t(fm1$alpha.mat)
fm1
#################################
## ordered factors (behaves like numerical variables):
data(soup, package="ordinal")
fm2 <- clm(SURENESS ~ 1, nominal=~PRODID + DAY, data=soup)
fm2$Theta
fm2$alpha.mat
prodid <- factor(soup$PRODID, ordered=TRUE)
fm2 <- clm(SURENESS ~ 1, nominal=~prodid + DAY, data=soup)
fm2$alpha.mat
fm2$Theta
fm2 <- clm(SURENESS ~ 1, nominal=~prodid, data=soup)
fm2$alpha.mat
fm2$Theta
#################################
## Aliased Coefficients:
##
## Example where the interaction in the nominal effects is aliased (by
## design). Here the two Theta matrices coincide. The alpha.mat
## matrices are similar except one has an extra row with NAs:
soup2 <- soup
levels(soup2$DAY)
levels(soup2$GENDER)
xx <- with(soup2, DAY == "2" & GENDER == "Female")
## Model with additive nominal effects:
fm8 <- clm(SURENESS ~ PRODID, nominal= ~ DAY + GENDER, data=soup2, subset=!xx)
fm8$alpha.mat
fm8$Theta
## Model with non-additive, but aliased nominal effects:
fm9 <- clm(SURENESS ~ PRODID, nominal= ~ DAY * GENDER, data=soup2, subset=!xx)
fm9$alpha.mat
fm9$Theta
stopEqual <- function(x, y, ca=FALSE)
stopifnot(isTRUE(all.equal(x, y, check.attributes=ca)))
stopEqual(fm8$alpha.mat, fm9$alpha.mat[1:3, ])
stopEqual(fm8$Theta, fm9$Theta)
stopEqual(logLik(fm8), logLik(fm9))
#################################
## Weights:
set.seed(12345)
wts <- runif(nrow(soup))
fm2 <- clm(SURENESS ~ 1, nominal=~SOUPTYPE + DAY, data=soup, weights=wts)
fm2$Theta
## Offset (correctly gives and error)
fm2 <- try(clm(SURENESS ~ 1, nominal=~SOUPTYPE + DAY + offset(wts),
data=soup), silent=TRUE)
stopifnot(inherits(fm2, "try-error"))
#################################
### Other (misc) examples:
fm2 <- clm(SURENESS ~ 1, nominal=~SOUPTYPE + DAY, data=soup)
fm2$Theta
fm2
fm2 <- clm(SURENESS ~ 1, nominal=~SOUPTYPE * DAY, data=soup)
fm2$Theta
fm2
fm2$alpha.mat
fm2 <- clm(SURENESS ~ 1, nominal=~SOUPTYPE * DAY, data=soup,
threshold="symmetric")
fm2$Theta
fm2$alpha.mat
#################################
### Check correctness of Theta matrix when intercept is removed in
### nominal formula:
### December 25th 2014, RHBC
fm1 <- clm(rating ~ temp, nominal=~contact-1, data=wine)
fm2 <- clm(rating ~ temp, nominal=~contact, data=wine)
stopifnot(isTRUE(all.equal(fm1$Theta, fm2$Theta)))
stopifnot(isTRUE(all.equal(fm1$logLik, fm2$logLik)))
wine2 <- wine
wine2$contact <- relevel(wine2$contact, "yes")
fm3 <- clm(rating ~ temp, nominal=~contact, data=wine2)
stopifnot(isTRUE(all.equal(coef(fm1, na.rm=TRUE), coef(fm3))))
#################################
|
