### NEGATIVE BINOMIAL
setwd("C:/Users/bha/Desktop/caroline/Data")
load('dat1.RData')
load('dat2.RData')
require(foreign)
## Loading required package: foreign
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.4
require(MASS)
## Loading required package: MASS
m1_nb <- glm.nb(sport ~ AGE + factor(Gend3) + factor(Eth7)+ factor(Disab3) +
factor(IMD10) + factor(BMIG) + factor(WorkStat10) +
factor(UrbRur2) + factor(CoastComm) + factor(Educ6), data = dat1)
summary(m1_nb)
##
## Call:
## glm.nb(formula = sport ~ AGE + factor(Gend3) + factor(Eth7) +
## factor(Disab3) + factor(IMD10) + factor(BMIG) + factor(WorkStat10) +
## factor(UrbRur2) + factor(CoastComm) + factor(Educ6), data = dat1,
## init.theta = 0.2880681274, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2170 -1.0553 -0.2933 0.1506 1.8356
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.095767 0.812849 9.960 < 2e-16 ***
## AGE -0.013357 0.005118 -2.610 0.00905 **
## factor(Gend3)2 -0.376546 0.147021 -2.561 0.01043 *
## factor(Eth7)2 -0.050454 0.432507 -0.117 0.90713
## factor(Disab3)2 0.259165 0.240546 1.077 0.28130
## factor(Disab3)3 0.107360 0.209178 0.513 0.60778
## factor(IMD10)2 -0.484812 0.390275 -1.242 0.21415
## factor(IMD10)3 -0.460689 0.404465 -1.139 0.25470
## factor(IMD10)4 -0.245640 0.385530 -0.637 0.52403
## factor(IMD10)5 -0.198316 0.407280 -0.487 0.62631
## factor(IMD10)6 -0.326299 0.383879 -0.850 0.39532
## factor(IMD10)7 -0.457439 0.383795 -1.192 0.23331
## factor(IMD10)8 -0.771968 0.346434 -2.228 0.02586 *
## factor(IMD10)9 -0.428588 0.355097 -1.207 0.22745
## factor(IMD10)10 -0.486668 0.365206 -1.333 0.18267
## factor(BMIG)2 -0.081670 0.581279 -0.141 0.88826
## factor(BMIG)3 -0.229846 0.585069 -0.393 0.69443
## factor(BMIG)4 -0.660538 0.595789 -1.109 0.26757
## factor(BMIG)5 -1.286469 0.727665 -1.768 0.07707 .
## factor(WorkStat10)2 -0.076774 0.219959 -0.349 0.72706
## factor(WorkStat10)5 -0.153533 0.177065 -0.867 0.38589
## factor(UrbRur2)2 0.050419 0.149835 0.336 0.73649
## factor(CoastComm)2 -0.205659 0.365948 -0.562 0.57412
## factor(Educ6)2 -0.281395 0.198065 -1.421 0.15540
## factor(Educ6)3 -0.320081 0.162835 -1.966 0.04934 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.2881) family taken to be 1)
##
## Null deviance: 919.24 on 726 degrees of freedom
## Residual deviance: 874.41 on 702 degrees of freedom
## AIC: 9588.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.2881
## Std. Err.: 0.0143
##
## 2 x log-likelihood: -9536.7250
### Logistric regression
m1_logit <- glm(bin ~ AGE + factor(Gend3) + factor(Eth7)+ factor(Disab3) +
factor(IMD10) + factor(BMIG) + factor(WorkStat10) +
factor(UrbRur2) + factor(CoastComm) + factor(Educ6), data = dat2, family = 'binomial')
summary(m1_logit)
##
## Call:
## glm(formula = bin ~ AGE + factor(Gend3) + factor(Eth7) + factor(Disab3) +
## factor(IMD10) + factor(BMIG) + factor(WorkStat10) + factor(UrbRur2) +
## factor(CoastComm) + factor(Educ6), family = "binomial", data = dat2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1693 -1.1505 0.6902 0.9539 1.7349
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.923427 0.993127 2.944 0.00324 **
## AGE -0.005810 0.005926 -0.980 0.32684
## factor(Gend3)2 -0.335387 0.172386 -1.946 0.05171 .
## factor(Eth7)2 -0.204983 0.510065 -0.402 0.68778
## factor(Disab3)2 0.408600 0.272595 1.499 0.13389
## factor(Disab3)3 0.560497 0.236215 2.373 0.01765 *
## factor(IMD10)2 -0.870451 0.507578 -1.715 0.08636 .
## factor(IMD10)3 -1.153150 0.516068 -2.234 0.02545 *
## factor(IMD10)4 -0.783376 0.501985 -1.561 0.11863
## factor(IMD10)5 -1.002176 0.524439 -1.911 0.05601 .
## factor(IMD10)6 -0.590789 0.504882 -1.170 0.24194
## factor(IMD10)7 -0.960778 0.498856 -1.926 0.05411 .
## factor(IMD10)8 -1.159322 0.459186 -2.525 0.01158 *
## factor(IMD10)9 -1.017552 0.468189 -2.173 0.02975 *
## factor(IMD10)10 -1.161093 0.476965 -2.434 0.01492 *
## factor(BMIG)2 -0.009844 0.685817 -0.014 0.98855
## factor(BMIG)3 -0.326054 0.688852 -0.473 0.63598
## factor(BMIG)4 -0.793658 0.700811 -1.132 0.25743
## factor(BMIG)5 -1.207148 0.856347 -1.410 0.15864
## factor(WorkStat10)2 -0.407564 0.254608 -1.601 0.10943
## factor(WorkStat10)5 -0.496034 0.205653 -2.412 0.01587 *
## factor(UrbRur2)2 0.202304 0.172904 1.170 0.24199
## factor(CoastComm)2 -0.695755 0.448390 -1.552 0.12074
## factor(Educ6)2 -0.327312 0.231228 -1.416 0.15691
## factor(Educ6)3 -0.560324 0.187377 -2.990 0.00279 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 972.78 on 726 degrees of freedom
## Residual deviance: 892.30 on 702 degrees of freedom
## AIC: 942.3
##
## Number of Fisher Scoring iterations: 4
### Multinomial baseline logit model (inactive 3 set as baseline)
library(nnet)
## Warning: package 'nnet' was built under R version 4.0.5
m1_multinom <- multinom(tri2 ~ AGE + factor(Gend3) + factor(Eth7)+ factor(Disab3) +
factor(IMD10) + factor(BMIG) + factor(WorkStat10) +
factor(UrbRur2) + factor(CoastComm) + factor(Educ6), data = dat2)
## # weights: 78 (50 variable)
## initial value 798.691134
## iter 10 value 630.123545
## iter 20 value 608.006303
## iter 30 value 605.798652
## iter 40 value 604.960082
## iter 50 value 604.853413
## iter 60 value 604.794542
## final value 604.793323
## converged
summary(m1_multinom)
## Call:
## multinom(formula = tri2 ~ AGE + factor(Gend3) + factor(Eth7) +
## factor(Disab3) + factor(IMD10) + factor(BMIG) + factor(WorkStat10) +
## factor(UrbRur2) + factor(CoastComm) + factor(Educ6), data = dat2)
##
## Coefficients:
## (Intercept) AGE factor(Gend3)2 factor(Eth7)2 factor(Disab3)2
## 1 2.451612 -0.01358481 -0.1359403 -0.7778084 0.6438725
## 2 -16.593287 -0.01928066 0.5142155 -1.8503381 0.7708757
## factor(Disab3)3 factor(IMD10)2 factor(IMD10)3 factor(IMD10)4 factor(IMD10)5
## 1 0.7294810 -0.4971983 -0.7862554 -0.2334186 -0.1917033
## 2 0.5742283 1.5055049 1.5288830 1.9935701 2.5584273
## factor(IMD10)6 factor(IMD10)7 factor(IMD10)8 factor(IMD10)9 factor(IMD10)10
## 1 0.1159537 -0.06567447 -0.7084951 -0.5462169 -0.9291285
## 2 2.3483587 2.69475414 1.7541741 1.7666894 0.8717315
## factor(BMIG)2 factor(BMIG)3 factor(BMIG)4 factor(BMIG)5 factor(WorkStat10)2
## 1 0.581015 0.1297968 -0.1526806 -0.457162 -0.42292843
## 2 13.655931 13.3084250 13.8009827 14.032336 -0.04850392
## factor(WorkStat10)5 factor(UrbRur2)2 factor(CoastComm)2 factor(Educ6)2
## 1 -0.49531311 0.03258696 -0.3067903 -0.5670284
## 2 -0.02086858 -0.46798269 1.6758582 -0.5814988
## factor(Educ6)3
## 1 -0.8379835
## 2 -0.7387098
##
## Std. Errors:
## (Intercept) AGE factor(Gend3)2 factor(Eth7)2 factor(Disab3)2
## 1 1.060495 0.007111597 0.2002827 0.5579242 0.3072612
## 2 1.344826 0.009456759 0.2823366 1.1263431 0.4237117
## factor(Disab3)3 factor(IMD10)2 factor(IMD10)3 factor(IMD10)4 factor(IMD10)5
## 1 0.2621320 0.5555947 0.5667544 0.556646 0.6096297
## 2 0.3735647 1.1956274 1.1848829 1.181006 1.2011819
## factor(IMD10)6 factor(IMD10)7 factor(IMD10)8 factor(IMD10)9 factor(IMD10)10
## 1 0.5736037 0.5738229 0.4982812 0.5102004 0.5116141
## 2 1.1810197 1.1676054 1.1288803 1.1356545 1.1671703
## factor(BMIG)2 factor(BMIG)3 factor(BMIG)4 factor(BMIG)5 factor(WorkStat10)2
## 1 0.7024913 0.7051726 0.7213461 0.9230374 0.3103067
## 2 0.4233403 0.4277507 0.4396651 0.6188948 0.4070757
## factor(WorkStat10)5 factor(UrbRur2)2 factor(CoastComm)2 factor(Educ6)2
## 1 0.2455830 0.2033663 0.4868654 0.2793210
## 2 0.3330804 0.2745518 1.1064480 0.3700823
## factor(Educ6)3
## 1 0.2216793
## 2 0.3042389
##
## Residual Deviance: 1209.587
## AIC: 1309.587
z <- summary(m1_multinom)$coefficients/summary(m1_multinom)$standard.errors
z
## (Intercept) AGE factor(Gend3)2 factor(Eth7)2 factor(Disab3)2
## 1 2.311762 -1.910234 -0.6787423 -1.394111 2.095522
## 2 -12.338615 -2.038823 1.8212852 -1.642784 1.819340
## factor(Disab3)3 factor(IMD10)2 factor(IMD10)3 factor(IMD10)4 factor(IMD10)5
## 1 2.782876 -0.8948939 -1.387295 -0.4193304 -0.3144586
## 2 1.537159 1.2591757 1.290324 1.6880274 2.1299249
## factor(IMD10)6 factor(IMD10)7 factor(IMD10)8 factor(IMD10)9 factor(IMD10)10
## 1 0.2021495 -0.1144508 -1.421878 -1.070593 -1.816073
## 2 1.9884162 2.3079321 1.553906 1.555657 0.746876
## factor(BMIG)2 factor(BMIG)3 factor(BMIG)4 factor(BMIG)5 factor(WorkStat10)2
## 1 0.8270778 0.1840638 -0.2116607 -0.4952801 -1.3629367
## 2 32.2575742 31.1125741 31.3897613 22.6732175 -0.1191521
## factor(WorkStat10)5 factor(UrbRur2)2 factor(CoastComm)2 factor(Educ6)2
## 1 -2.01688678 0.1602378 -0.6301338 -2.030024
## 2 -0.06265327 -1.7045331 1.5146290 -1.571269
## factor(Educ6)3
## 1 -3.780162
## 2 -2.428059
# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
## (Intercept) AGE factor(Gend3)2 factor(Eth7)2 factor(Disab3)2
## 1 0.02079078 0.05610313 0.49730116 0.1632840 0.03612462
## 2 0.00000000 0.04146767 0.06856352 0.1004277 0.06885958
## factor(Disab3)3 factor(IMD10)2 factor(IMD10)3 factor(IMD10)4 factor(IMD10)5
## 1 0.005387933 0.3708438 0.1653519 0.67497471 0.75317278
## 2 0.124254411 0.2079669 0.1969381 0.09140596 0.03317781
## factor(IMD10)6 factor(IMD10)7 factor(IMD10)8 factor(IMD10)9 factor(IMD10)10
## 1 0.83979981 0.90888047 0.1550617 0.2843525 0.06935919
## 2 0.04676567 0.02100291 0.1202068 0.1197896 0.45513843
## factor(BMIG)2 factor(BMIG)3 factor(BMIG)4 factor(BMIG)5 factor(WorkStat10)2
## 1 0.408193 0.8539634 0.8323718 0.6204024 0.1729025
## 2 0.000000 0.0000000 0.0000000 0.0000000 0.9051549
## factor(WorkStat10)5 factor(UrbRur2)2 factor(CoastComm)2 factor(Educ6)2
## 1 0.04370733 0.87269377 0.5286071 0.04235408
## 2 0.95004261 0.08828155 0.1298664 0.11612019
## factor(Educ6)3
## 1 0.0001567266
## 2 0.0151798894