#Bryce Sellner Quiz 2
library(faraway)
data(africa)
?africa
newafrica <- na.omit(africa)
attach(newafrica)
#checking to see if poisson regression is called for
hist(miltcoup)
#the histogram is right skewed and looks like that of a poisson pmf. This means using poisson regression is appropriate
lmc.o <- glm(miltcoup ~ oligarchy, family = "poisson")
lmc.pol <- glm(miltcoup ~ pollib, family = "poisson")
lmc.pa <- glm(miltcoup ~ parties , family = "poisson")
lmc.pc <- glm(miltcoup ~ pctvote , family = "poisson")
lmc.pop <- glm(miltcoup ~ popn , family = "poisson")
lmc.s <- glm(miltcoup ~ size, family = "poisson")
lmc.nume <- glm(miltcoup ~ numelec, family = "poisson")
lmc.numr <- glm(miltcoup ~ numregim, family = "poisson")
AIC(lmc.o)
## [1] 111.7554
AIC(lmc.pol)
## [1] 127.8726
AIC(lmc.pa )
## [1] 129.1635
AIC(lmc.pc )
## [1] 134.748
AIC(lmc.pop)
## [1] 128.5262
AIC(lmc.s )
## [1] 133.5853
AIC(lmc.nume )
## [1] 134.691
AIC(lmc.numr )
## [1] 130.1283
#Checking the AICs, the number of years that a country was ruled by a military oligarchy is the best predictor of number of military coups since it was the smallest by a decent margin
coef(lmc.o)
## (Intercept) oligarchy
## -0.2544511 0.0998016
exp(0.0998016) #Taking the exponent to find out how number of years in a oligarchy affects number of military coups
## [1] 1.104952
#An increase of one extra year under a military coup leads to a multiplicative increase of 1.105 in the number of military coups a country has
#Start of drop in deviance tests, making a model for each other variable + oligarchy. I will test them using anove and a significance level of .05
lmc.o.pol <- glm(miltcoup ~ oligarchy + pollib, family = "poisson")
lmc.o.pa <- glm(miltcoup ~ oligarchy + parties, family = "poisson")
lmc.o.pc <- glm(miltcoup ~ oligarchy + pctvote, family = "poisson")
lmc.o.pop <- glm(miltcoup ~ oligarchy + popn, family = "poisson")
lmc.o.s <- glm(miltcoup ~ oligarchy + size, family = "poisson")
lmc.o.nume <- glm(miltcoup ~ oligarchy + numelec, family = "poisson")
lmc.o.numr <- glm(miltcoup ~ oligarchy + numregim, family = "poisson")
anova(lmc.o, lmc.o.pol, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + pollib
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 38.204 1 4.7429 0.02942 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmc.o, lmc.o.pa, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + parties
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 40.025 1 2.9213 0.08742 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmc.o, lmc.o.pc, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + pctvote
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 42.914 1 0.032627 0.8567
anova(lmc.o, lmc.o.pop, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + popn
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 42.705 1 0.24124 0.6233
anova(lmc.o, lmc.o.s, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + size
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 42.805 1 0.14161 0.7067
anova(lmc.o, lmc.o.nume, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + numelec
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 39.719 1 3.2278 0.0724 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmc.o, lmc.o.numr, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy
## Model 2: miltcoup ~ oligarchy + numregim
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 34 42.947
## 2 33 42.760 1 0.18679 0.6656
#Anova is used to perform the actual drop in deviance test
#Only one of the new models has a pvalue under .05, that being the model with political liberalization added to it
#Next is another step of of drop in deviance tests, this time adding each variable to both oligarchy and pollib
lmc.o.pol.pa <- glm(miltcoup ~ oligarchy + pollib + parties, family = "poisson")
lmc.o.pol.pc <- glm(miltcoup ~ oligarchy + pollib + pctvote, family = "poisson")
lmc.o.pol.pop <- glm(miltcoup ~ oligarchy + pollib + popn , family = "poisson")
lmc.o.pol.s <- glm(miltcoup ~ oligarchy + pollib + size, family = "poisson")
lmc.o.pol.nume <- glm(miltcoup ~ oligarchy + pollib + numelec , family = "poisson")
lmc.o.pol.numr <- glm(miltcoup ~ oligarchy + pollib + numregim , family = "poisson")
anova(lmc.o.pol, lmc.o.pol.pa, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + parties
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 32.856 1 5.348 0.02075 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmc.o.pol, lmc.o.pol.pc, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + pctvote
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 37.599 1 0.60461 0.4368
anova(lmc.o.pol, lmc.o.pol.pop, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + popn
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 38.198 1 0.0055912 0.9404
anova(lmc.o.pol, lmc.o.pol.s, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + size
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 37.576 1 0.62727 0.4284
anova(lmc.o.pol, lmc.o.pol.nume, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + numelec
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 37.147 1 1.0567 0.304
anova(lmc.o.pol, lmc.o.pol.numr, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib
## Model 2: miltcoup ~ oligarchy + pollib + numregim
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 33 38.204
## 2 32 37.749 1 0.45472 0.5001
#Only one of the new models has a pvalue under .05, that being the model with number of legal political parties added to it
#Next is another step of of drop in deviance tests, this time adding each variable to oligarchy, pollib, and parties
lmc.o.pol.pa.pc <- glm(miltcoup ~ oligarchy + pollib + parties + pctvote, family = "poisson")
lmc.o.pol.pa.pop <- glm(miltcoup ~ oligarchy + pollib + parties + popn, family = "poisson")
lmc.o.pol.pa.s <- glm(miltcoup ~ oligarchy + pollib + parties + size, family = "poisson")
lmc.o.pol.pa.nume <- glm(miltcoup ~ oligarchy + pollib + parties + numelec, family = "poisson")
lmc.o.pol.pa.numr <- glm(miltcoup ~ oligarchy + pollib + parties + numregim, family = "poisson")
anova(lmc.o.pol.pa, lmc.o.pol.pa.pc, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib + parties
## Model 2: miltcoup ~ oligarchy + pollib + parties + pctvote
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 32.856
## 2 31 31.081 1 1.7752 0.1827
anova(lmc.o.pol.pa, lmc.o.pol.pa.pop, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib + parties
## Model 2: miltcoup ~ oligarchy + pollib + parties + popn
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 32.856
## 2 31 32.247 1 0.60899 0.4352
anova(lmc.o.pol.pa, lmc.o.pol.pa.s, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib + parties
## Model 2: miltcoup ~ oligarchy + pollib + parties + size
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 32.856
## 2 31 32.539 1 0.31715 0.5733
anova(lmc.o.pol.pa, lmc.o.pol.pa.nume, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib + parties
## Model 2: miltcoup ~ oligarchy + pollib + parties + numelec
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 32.856
## 2 31 32.689 1 0.16705 0.6828
anova(lmc.o.pol.pa, lmc.o.pol.pa.numr, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: miltcoup ~ oligarchy + pollib + parties
## Model 2: miltcoup ~ oligarchy + pollib + parties + numregim
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32 32.856
## 2 31 32.692 1 0.16329 0.6861
#None of the models in this run had a significant p-value, therefore none of the other variables should be added. This leaves lmc.o.pol.pa as the Final model
coef(lmc.o.pol.pa)
## (Intercept) oligarchy pollib parties
## 0.25137715 0.09262231 -0.57410330 0.02205864
exp(0.09262231)
## [1] 1.097047
exp(-0.57410330)
## [1] 0.5632097
exp(0.02205864) #Again taking exponents of parameters to see their effect on the response
## [1] 1.022304
#Regarding the final model, an increase of one extra year under a military coup leads to a multiplicative increase of 1.097 in the number of military coups a country has, an increase of one extra legal political party leads to a multiplicative increase of 1.022 in the number of military coups a country has, and each category of higher political liberalization leads to a multiplicative decrease of .563 in the number of military coups
#Now, I will check for a 95% confidence interval of the effect the number of political parties has on miliary coups
exp(confint(lmc.o.pol.pa, parm="parties", level = .95))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## 1.003563 1.039670
#Using the confint function to find the 95% confidence interval of the parameter "parties" and then taking the exponents of them, I am 95% confident that an increase in one legal political party will lead to a multiplicative increase of 1.0036 to 1.0397 in the number of military coups
#Finally, Ill test for overdispersion to see how well the model measures the variance in the response
mcdp<-sum(residuals(lmc.o.pol.pa,type="pearson")^2/lmc.o.pol.pa$df.res)
mcdp
## [1] 0.9841373
sqrt(mcdp)
## [1] 0.992037
#taking the squareroot of the dispersion parameter yields .992. This means the standard errors need to be multplied by .992. This is a good sign since they barely have to be changed and means the model did an accurate job