-
Notifications
You must be signed in to change notification settings - Fork 15
Expand file tree
/
Copy pathCh1-Introduction.R
More file actions
395 lines (339 loc) · 12.2 KB
/
Copy pathCh1-Introduction.R
File metadata and controls
395 lines (339 loc) · 12.2 KB
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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
#-----------------------------------------------------------------
# Chapter 1 - Introduction
# Computational aspects of psychometric methods. With R.
# P. Martinkova & A. Hladka
#-----------------------------------------------------------------
#-----------------------------------------------------------------
# Plot settings
#-----------------------------------------------------------------
theme_fig <- function(base_size = 17, base_family = "") {
theme_bw(base_size = base_size, base_family = base_family) +
theme(
legend.key = element_rect(fill = "white", colour = NA),
axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank(),
legend.background = element_blank()
)
}
#-----------------------------------------------------------------
# 1.5 Exploring measurement data
#-----------------------------------------------------------------
#-----------------------------------------------------------------
# 1.5.1 Item scores
#-----------------------------------------------------------------
#-----------------------------------------------------------------
# Nominal items
#-----------------------------------------------------------------
#--------------
# loading data
data(HCItest, package = "ShinyItemAnalysis")
# view data head
head(HCItest, n = 2)
## Item 1 Item 2 Item 3 Item 4 Item 5 Item 6 Item 7 Item 8 Item 9 ...
## 1 D B A D B B B C D ...
## 2 D B A D B C B C D ...
#--------------
#--------------
table(HCItest$"Item 1")
## A B C D
## 27 59 110 455
proportions(table(HCItest$"Item 1"))
prop.table(table(HCItest$"Item 1"))
## A B C D
## 0.0415 0.0906 0.1690 0.6989
#--------------
#--------------
data("HCIkey", package = "ShinyItemAnalysis")
unlist(HCIkey)
## key1 key2 key3 key4 key5 key6 key7 key8 key9 key10 key11
## D B A D B C C C D A A
## key12 key13 key14 key15 key16 key17 key18 key19 key20
## D A A C A C C C D
## Levels: A B C D
#--------------
#--------------
HCIscored <- as.data.frame(mirt::key2binary(HCItest[, 1:20], HCIkey))
head(HCIscored, n = 2)
## Item 1 Item 2 Item 3 Item 4 Item 5 Item 6 Item 7 Item 8 Item 9 ...
## 1 1 1 1 1 1 0 0 1 1 ...
## 2 1 1 1 1 1 1 0 1 1 ...
#--------------
#--------------
data(HCI, package = "ShinyItemAnalysis")
head(HCI, n = 2)
## Item 1 Item 2 Item 3 Item 4 Item 5 Item 6 Item 7 Item 8 Item 9 ...
## 1 1 1 1 1 1 0 0 1 1 ...
## 2 1 1 1 1 1 1 0 1 1 ...
#--------------
#-----------------------------------------------------------------
# Binary items
#-----------------------------------------------------------------
#--------------
table(HCI$"Item 1")
## 0 1
## 196 455
#--------------
#--------------
proportions(table(HCI$"Item 1"))
prop.table(table(HCI$"Item 1"))
## 0 1
## 0.3011 0.6989
#--------------
#--------------
summary(HCI$"Item 1")
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.6989 1.0000 1.0000
#--------------
#--------------
# summary statistics for all variables (not shown in the book)
summary(HCI)
## Item1 Item2 Item3 Item4
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:0.000
## Median :1.0000 Median :1.0000 Median :1.0000 Median :0.000
## Mean :0.6989 Mean :0.7527 Mean :0.8479 Mean :0.404
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.000
## ...
#--------------
#-----------------------------------------------------------------
# Ordinal items
#-----------------------------------------------------------------
#--------------
data("BFI2", package = "ShinyItemAnalysis")
head(BFI2, n = 2)
## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12 i13 i14 i15 i16 i17 i18 ...
## 1 5 5 2 3 5 4 4 3 4 5 5 4 3 4 2 3 5 3 ...
## 2 4 5 4 3 3 3 5 3 2 4 4 2 4 2 3 4 5 4 ...
#--------------
#--------------
table(BFI2$i1)
## 1 2 3 4 5
## 50 264 321 758 340
#--------------
#--------------
BFI2binary <- 1 * as.data.frame(BFI2[, 1:60] >= 3)
head(BFI2binary, n = 2)
## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12 i13 i14 i15 i16 i17 i18 ...
## 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 ...
## 2 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 ...
#--------------
#-----------------------------------------------------------------
# Continuous items
#-----------------------------------------------------------------
#--------------
data(EPIA, package = "ShinyItemAnalysis")
head(EPIA, n = 2)
## Item 1 Item 2 Item 3 Item 4 Item 5 score
## 1 96 36 80 78 79 369
## 2 42 2 1 1 1 47
#--------------
#--------------
summary(EPIA$"Item 1")
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0000 40.0000 63.0000 62.8296 89.0000 111.0000
#--------------
#-----------------------------------------------------------------
# 1.5.2 Test scores
#-----------------------------------------------------------------
#--------------
HCI$total
# [1] 16 19 17 20 19 20 20 14 18 17 17 16 15 12 17 ...
rowSums(HCI[, 1:20])
# [1] 16 19 17 20 19 20 20 14 18 17 17 16 15 12 17 ...
apply(X = HCI[, 1:20], MARGIN = 1, FUN = sum)
# [1] 16 19 17 20 19 20 20 14 18 17 17 16 15 12 17 ...
#--------------
#--------------
# summary of total score
summary(HCI$total)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.0000 10.0000 12.0000 12.2120 15.0000 20.0000
#--------------
#--------------
c(Min = min(HCI$total), Max = max(HCI$total), Mean = mean(HCI$total),
Med = median(HCI$total), Var = var(HCI$total), SD = sd(HCI$total),
Skew = moments::skewness(HCI$total), Kurt = moments::kurtosis(HCI$total))
## Min Max Mean Med Var SD Skew Kurt
## 3.0000 20.0000 12.2120 12.0000 13.2473 3.6397 -0.1982 2.3474
#--------------
#--------------
psych::describe(HCI$total, type = 1)
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 651 12.21 3.64 12 12.32 4.45 3 20 17 -0.2 -0.65 0.14
#--------------
#--------------
# obtaining the same value for kurtosis, not shown in the book
moments::kurtosis(HCI$total) - 3
## [1] -0.6526
#--------------
#-----------------------------------------------------------------
# Standardizations and other transformations
#-----------------------------------------------------------------
#--------------
# Z-score
(zscore <- as.vector(scale(HCI$total)))
## [1] 1.0408 1.8650 1.3155 2.1398 1.8650 2.1398 2.1398 ...
# T-score
(tscore <- 10 * zscore + 50)
## [1] 60.4075 68.6500 63.1550 71.3975 68.6500 71.3975 71.3975 ...
# success rate
(success_rate <- 100 * (HCI$total / max(HCI$total)))
## [1] 80 95 85 100 95 100 100 ...
#--------------
#--------------
plot(ecdf(HCI$total), xlab = "HCI score", ylab = "Percentile")
ecdf(HCI$total)(HCI$total)
## [1] 0.8725 0.9923 0.9401 1.0000 0.9923 1.0000 1.0000 ...
(percentiles <- round(100 * ecdf(HCI$total)(HCI$total))) # percentiles
## [1] 87 99 94 100 99 100 100 ...
#--------------
#--------------
# Summarizing the score, its standardizations and transformations
# (code not displayed in the book)
head(data.frame(score = HCI$total, zscore, tscore, percentiles, success_rate), n = 4)
## score zscore tscore percentiles success_rate
## 1 16 1.0408 60.4075 87 80
## 2 19 1.8650 68.6500 99 95
## 3 17 1.3155 63.1550 94 85
## 4 20 2.1398 71.3975 100 100
#--------------
#-----------------------------------------------------------------
# 1.5.3 Covariates
#-----------------------------------------------------------------
#--------------
table(HCI$"Item 1", HCI$major)
## 0 1
## 0 97 99
## 1 168 287
#--------------
#--------------
proportions(table(HCI$"Item 1", HCI$major), margin = 2)
prop.table(table(HCI$"Item 1", HCI$major), margin = 2)
## 0 1
## 0 0.3660 0.2565
## 1 0.6340 0.7435
#--------------
#--------------
table(HCItest$gender, HCItest$"Item 1")
## A B C D
## 0 16 34 65 290
## 1 11 25 45 165
proportions(table(HCItest$gender, HCItest$"Item 1"), margin = 1)
prop.table(table(HCItest$gender, HCItest$"Item 1"), margin = 1)
## A B C D
## 0 0.0395 0.0840 0.1605 0.7160
## 1 0.0447 0.1016 0.1829 0.6707
by(HCI$"Item 1", HCI$gender, FUN = summary)
## HCI$gender: 0
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 0.716 1.000 1.000
## ----------------------------------------------------
## HCI$gender: 1
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.6707 1.0000 1.0000
#--------------
#-----------------------------------------------------------------
# 1.6 Modelling measurement data
#-----------------------------------------------------------------
#-----------------------------------------------------------------
# 1.6.1 Discrete random variables
#-----------------------------------------------------------------
#--------------
# Bernoulli distribution
dbinom(x = 1, size = 1, prob = 0.3)
## [1] 0.3
dbinom(x = 0, size = 1, prob = 0.3)
## [1] 0.7
#--------------
#--------------
# simulating a Bernoulli trial
set.seed(42)
(item1 <- rbinom(n = 100, size = 1, prob = 0.3))
## [1] 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 1 0 0 0 0 0 1 0 1 1 1 1 0 0 0 0 0
## [37] 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0
## [73] 1 0 1 1 0 0 0 1 0 0 0 0 0 1 0 1 0 0 1 0 0 0 1 0 1 1 0 0
# sample mean and variance
mean(item1)
## [1] 0.3400
var(item1)
## [1] 0.2267
#--------------
#--------------
# probability of gaining total score of exact value
dbinom(5, size = 10, prob = 0.5)
## [1] 0.2461
dbinom(10, size = 10, prob = 0.5)
## [1] 0.0010
#--------------
#--------------
# probability of gaining total score at least of given value or in range
pbinom(5, size = 10, prob = 0.5)
## [1] 0.6230
pbinom(10, size = 10, prob = 0.5)
## [1] 1.0000
pbinom(8, size = 10, prob = 0.5) - pbinom(2, size = 10, prob = 0.5)
## [1] 0.9346
#--------------
#--------------
# simulating from a binomial distribution
score <- rbinom(n = 100, size = 20, prob = 0.7)
mean(score)
## [1] 13.9600
20 * 0.7
## [1] 14.0000
sd(score)
## [1] 1.9065
sqrt(20 * 0.7 * (1 - 0.7))
## [1] 2.0494
#--------------
#-----------------------------------------------------------------
# 1.6.2 Continuous random variables
#-----------------------------------------------------------------
#--------------
dnorm(x = 5, mean = 10 * 0.5, sd = sqrt(10 * 0.5 * 0.5))
## [1] 0.2523
dnorm(x = 10, mean = 10 * 0.5, sd = sqrt(10 * 0.5 * 0.5))
## [1] 0.0017
pnorm(q = 5, mean = 10 * 0.5, sd = sqrt(10 * 0.5 * 0.5))
## [1] 0.5000
pnorm(q = 10, mean = 10 * 0.5, sd = sqrt(10 * 0.5 * 0.5))
## [1] 0.9992
#--------------
#--------------
qnorm(0.975)
## [1] 1.9600
qnorm(0.025)
## [1] -1.9600
#--------------
#-------------
# Histogram with an estimated normal density
library(ggplot2)
ggplot(data = HCI, aes(total)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 1,
col = "black", fill = "gold") +
stat_function(fun = dnorm, colour = "red", linewidth = 0.8,
args = list(mean = mean(HCI$total), sd = sd(HCI$total))) +
xlab("Total score") + ylab("Density") + theme_fig()
#--------------
#--------------
# QQ plot in base
qqnorm(HCI$total)
qqline(HCI$total)
#--------------
#--------------
# QQ plot in ggplot
ggplot(HCI, aes(sample = total)) +
stat_qq(size = 2, shape = 1) +
stat_qq_line() +
theme_fig() +
ylab("Sample Quantiles") +
xlab("Theoretical Quantiles")
#--------------
#-----------------------------------------------------------------
# 1.7 ShinyItemAnalysis interactive application
#-----------------------------------------------------------------
ShinyItemAnalysis::run_app()