1102class1
2 input data
## Score Nomination
## 1 1.56 0
## 2 1.56 0
## 3 1.11 0
## 4 1.56 0
## 5 1.22 4
## 6 1.33 0## 'data.frame': 291 obs. of 2 variables:
## $ Score : num 1.56 1.56 1.11 1.56 1.22 1.33 1.11 1 1.11 1.22 ...
## $ Nomination: int 0 0 0 0 4 0 0 0 0 19 ...3 data management
## [1] 0.5635739## Score Nomination
## 1.658110 2.487973## Score Nomination
## 0.4842037 41.26451004 histogram of nomination
ggplot(dta, aes(Nomination)) +
stat_bin(binwidth = bd) +
labs(x = "Number of nominations", y = "Count")
# poisson fit
ggplot(dta, aes(Score, Nomination)) +
geom_jitter(alpha = 0.2) +
stat_smooth(method = "glm", method.args = list(family = poisson)) +
labs(y = "Number of peer nomination", x = "Bully score")## `geom_smooth()` using formula 'y ~ x'
從圖中可看到，被同儕視為霸凌者的資料分布，多數的人都不是霸凌者(0的數量很高)，但這些0裡面有依些是真的0(真的不是霸凌者)以及是假的0(是潛在的霸凌者)
用poisson去fit bully score 和 Number of peer nomination的結果圖顯示，bully score高的人也比較會被同儕認為是霸凌者
5 Model
##
## Call:
## glm(formula = Nomination ~ Score, family = poisson, data = dta)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.4995 -1.8246 -1.5952 -0.1552 11.6806
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.66308 0.08871 -7.475 7.75e-14 ***
## Score 0.81434 0.03504 23.239 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2205.2 on 290 degrees of freedom
## Residual deviance: 1774.6 on 289 degrees of freedom
## AIC: 2160.3
##
## Number of Fisher Scoring iterations: 6##
## Call:
## zeroinfl(formula = Nomination ~ Score | Score, data = dta)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -2.37068 -0.68873 -0.59685 -0.05357 9.87491
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.49604 0.09987 4.967 6.81e-07 ***
## Score 0.59609 0.03942 15.120 < 2e-16 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4967 0.3463 4.322 1.54e-05 ***
## Score -0.7716 0.1967 -3.924 8.72e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 8
## Log-likelihood: -780.6 on 4 Df# zero-inflated negative binomial
summary(m2 <- zeroinfl(Nomination ~ Score | Score, data = dta,
dist = "negbin"))##
## Call:
## zeroinfl(formula = Nomination ~ Score | Score, data = dta, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.58338 -0.48698 -0.38705 0.02297 6.96982
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6737 0.4354 -1.547 0.122
## Score 0.8819 0.2052 4.297 1.73e-05 ***
## Log(theta) -1.0658 0.1792 -5.947 2.73e-09 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.954 2.484 1.189 0.234
## Score -3.419 2.260 -1.513 0.130
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 0.3445
## Number of iterations in BFGS optimization: 21
## Log-likelihood: -494.7 on 5 DfIn m0 poisson(no Zero-inflation),分數增加一單位，被同儕認為是霸凌者的次數增加[exp(0.814)-1]x100%。
In m1 poisson,分數增加一單位，被同儕認為是霸凌者的次數增加exp(0.596)。在沒有被同儕認為是霸凌者中，分數增加一單位，被同儕認為是霸凌者的次數減少exp(0.7716)。
In m2 negbin,分數增加一單位，被同儕認為是霸凌者的次數增加exp(0.8819)。在沒有被同儕認為是霸凌者中，分數增加一單位，被同儕認為是霸凌者的次數減少exp(3.419)。
6 Model comparison
## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -4.891816 model2 > model1 4.9955e-07
## AIC-corrected -4.858930 model2 > model1 5.9011e-07
## BIC-corrected -4.798531 model2 > model1 7.9917e-07## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -4.242754 model2 > model1 1.104e-05
## AIC-corrected -4.242754 model2 > model1 1.104e-05
## BIC-corrected -4.242754 model2 > model1 1.104e-05兩兩模型比較，皆有顯著。比較兩個模型比較結果，在H_A欄(model2 > model1)，所以看起來是用negative binomial的model fit 比較好
7 fortify data with fitted values and residuals
dta_m <- data.frame(dta, fit1 = predict(m1, type = "response"),
fit2 = predict(m2, type = "response"),
r1 = resid(m1, type="pearson"),
r2 = resid(m2, type="pearson"))
# residual plots
ggplot(dta_m, aes(fit1, r1)) +
geom_point(pch=4) +
geom_point(aes(fit2, r2), pch=1) +
geom_hline(yintercept = 0) +
labs(x = "Fitted values", y = "Residuals",
title = "Poisson vs. Negative Binomial")
用圖比較poisson和negbin兩個 fitted values and residuals的關係，看起來的確是negbin的值相對離residual=0比較近
# fortify data with model fits and CIs
dta_m2 <- data.frame(dta, yhat = predict(m2))
# model fit over observations
ggplot(dta_m2, aes(Score, Nomination)) +
geom_point(pch = 20, alpha = .2) +
geom_line(aes(Score, yhat), col = "magenta", lwd=rel(1)) +
stat_smooth(method = "glm", method.args = list(family = poisson)) +
labs(x = "Bully Score", y = "Predicted number of nomination")## `geom_smooth()` using formula 'y ~ x'
桃紅色線代表negbin的fitted value curve，藍色線為poisson的fitted value curve
兩個curve顯示，bully score越高被同儕認為是霸凌者也越高，且用negbin估計的比poisson高