-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathWorkshop_code_along.R
173 lines (145 loc) · 4.4 KB
/
Workshop_code_along.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
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
install.packages("singleRcapture")
remotes::install_github("https://github.com/ncn-foreigners/singleRcaptureExtra")
library(singleRcapture)
library(singleRcaptureExtra)
df_drunk <- read.csv("~/Desktop/Workshop-on-Survey-Statistics-2024-singleRcapture/df_drunk")
head(df_drunk)
tab <- table(df_drunk$citizenship)
nmm <- names(tab[tab > 200])
df_drunk$citizenship <- ifelse(
df_drunk$citizenship %in% nmm,
df_drunk$citizenship,
"OTHER"
)
#library(countreg)
model <- estimatePopsize(
formula = counts ~ gender + citizenship,
data = df_drunk,
model = singleRcapture::ztpoisson()
)
summary(model)
model_2 <- estimatePopsize(
formula = counts ~ gender + citizenship,
data = df_drunk,
model = ztoipoisson(omegaLink = "cloglog"),
controlModel = controlModel(
omegaFormula = ~ gender + citizenship
),
controlMethod = controlMethod(
verbose = 5
)
)
summary(model_2)
model_2_reduced <- estimatePopsize(
formula = counts ~ 1,
data = df_drunk,
model = ztoipoisson(omegaLink = "cloglog"),
controlModel = controlModel(
omegaFormula = ~ gender + poly(age, degree = 3)
),
controlMethod = controlMethod(
verbose = 5
)
)
summary(model_2_reduced)
AIC(model_2)
AIC(model_2_reduced)
library(lmtest)
lrtest(model_2, model_2_reduced)
lrtest(model_2_reduced, model)
model_3 <- estimatePopsize(
formula = counts ~ 1,
data = df_drunk,
model = ztoinegbin(
omegaLink = "cloglog"
),
controlModel = controlModel(
omegaFormula = ~ gender + poly(age, degree = 3)
),
controlMethod = controlMethod(
verbose = 5, stepsize = .25
)
)
summary(model_3)
summary(marginalFreq(model_2_reduced), dropl5 = "group", df = 1)
summary(marginalFreq(model_3), dropl5 = "group", df = 1)
lrtest(model_2_reduced, model_3)
library(VGAM)
library(VGAMdata)
additive_model <- vgam(
formula = counts ~ gender + sm.os(age, df = 3),
data = df_drunk,
family = VGAMdata::oipospoisson(lpstr1 = "clogloglink"),
trace = TRUE
)
model_4 <- estimatePopsize(additive_model)
plot(additive_model, which.term = 2)
summary(model_4)
AIC(model_3)
AIC(model_4)
# model_5 <- estimatePopsize(
# formula = counts ~ gender + age + citizenship,
# data = df_drunk,
# model = singleRcapture::ztpoisson,
# popVar = "bootstrap",
# controlPopVar = controlPopVar(
# bootType = "semiparametric",
# confType = "basic",
# alpha = .001,
# traceBootstrapSize = TRUE,
# bootstrapVisualTrace = TRUE
# )
# )
model_5 <- estimatePopsize(
formula = counts ~ 1,
data = df_drunk,
model = Hurdleztgeom(),
controlModel = controlModel(
piFormula = ~ gender + poly(age, degree = 3)
)
)
model_6 <- estimatePopsize(
formula = counts ~ 1,
data = df_drunk,
model = Hurdleztnegbin(),
controlModel = controlModel(
piFormula = ~ gender + poly(age, degree = 3)
),
controlMethod = controlMethod(
verbose = 4,
stepsize = .35,
momentumFactor = 1
)
)
?controlMethod
library(ggplot2)
fitt_data <- data.frame(
"point" = c(model$populationSize$pointEstimate,
model_2_reduced$populationSize$pointEstimate,
model_3$populationSize$pointEstimate,
model_4$populationSize$pointEstimate,
model_5$populationSize$pointEstimate,
model_6$populationSize$pointEstimate),
"lower" = c(model$populationSize$confidenceInterval[1, 1],
model_2_reduced$populationSize$confidenceInterval[1, 1],
model_3$populationSize$confidenceInterval[1, 1],
model_4$populationSize$confidenceInterval[1, 1],
model_5$populationSize$confidenceInterval[1, 1],
model_6$populationSize$confidenceInterval[1, 1]),
"upper" = c(model$populationSize$confidenceInterval[1, 2],
model_2_reduced$populationSize$confidenceInterval[1, 2],
model_3$populationSize$confidenceInterval[1, 2],
model_4$populationSize$confidenceInterval[1, 2],
model_5$populationSize$confidenceInterval[1, 2],
model_6$populationSize$confidenceInterval[1, 2]),
"name" = c("ztpoisson", "ztoipoisson", "ztoinegbin", "additive",
"Hurdleztgeom", "Hurdleztnegbin")
)
library(ggplot2)
ggplot(fitt_data) +
geom_point(aes(y = name, x = point)) +
geom_errorbar(aes(y = name, x = point, xmin = lower, xmax = upper))
# dpop <- dfpopsize(model_3, cores = 5)
# plot(model_3, plotType = "dfpopContr", )
plot(model_3, plotType = "scaleLoc")
?singleRcapture:::plot.singleRStaticCountData