Homework 9
Abby_33
April 25, 2015
library("Zelig")## Loading required package: MASS
## Loading required package: boot## ##
## ## Zelig (Version 3.5.3, built: 2011-11-29)
## ## Please refer to http://gking.harvard.edu/zelig for full
## ## documentation or help.zelig() for help with commands and
## ## models supported by Zelig.
## ##
##
## ## Zelig project citations:
## ## Kosuke Imai, Gary King, and Olivia Lau. (2009).
## ## ``Zelig: Everyone's Statistical Software,''
## ## http://gking.harvard.edu/zelig.
## ## and
## ## Kosuke Imai, Gary King, and Olivia Lau. (2008).
## ## ``Toward A Common Framework for Statistical Analysis
## ## and Development,'' Journal of Computational and
## ## Graphical Statistics, Vol. 17, No. 4 (December)
## ## pp. 892-913.
##
## ## To cite individual Zelig models, please use the citation format printed with
## ## each model run and in the documentation.
## ##library("DescTools")
library("dplyr")## Warning: package 'dplyr' was built under R version 3.1.3##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:MASS':
##
## select
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("stargazer")##
## Please cite as:
##
## Hlavac, Marek (2014). stargazer: LaTeX code and ASCII text for well-formatted regression and summary statistics tables.
## R package version 5.1. http://CRAN.R-project.org/package=stargazerlibrary("readstata13")## Warning: package 'readstata13' was built under R version 3.1.3library("foreign")
library("car")##
## Attaching package: 'car'
##
## The following object is masked from 'package:DescTools':
##
## Recode
##
## The following object is masked from 'package:boot':
##
## logitlibrary("npmlreg")## Warning: package 'npmlreg' was built under R version 3.1.3library("texreg")## Warning: package 'texreg' was built under R version 3.1.3## Version: 1.34
## Date: 2014-10-31
## Author: Philip Leifeld (University of Konstanz)
##
## Please cite the JSS article in your publications -- see citation("texreg").library("lme4")## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following object is masked from 'package:npmlreg':
##
## expand
##
## Loading required package: Rcpplibrary("magrittr")## Warning: package 'magrittr' was built under R version 3.1.3##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:texreg':
##
## extractdata(irlsuicide)str(irlsuicide)## 'data.frame': 104 obs. of 8 variables:
## $ Region : Factor w/ 13 levels "Cork ","Dublin ",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ ID : Factor w/ 13 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ pop : int 31923 31907 8517 8127 13527 12965 11958 8299 117575 114626 ...
## $ death : int 6 52 8 37 17 33 14 22 33 114 ...
## $ sex : Factor w/ 2 levels "0","1": 1 2 1 2 1 2 1 2 1 2 ...
## $ age : Factor w/ 4 levels "1","2","3","4": 1 1 2 2 3 3 4 4 1 1 ...
## $ smr : num 97.7 144 181.7 193.4 182.4 ...
## $ expected: num 6.14 36.11 4.4 19.13 9.32 ...summary(irlsuicide)## Region ID pop death sex
## Cork : 8 1 : 8 Min. : 2434 Min. : 1.00 0:52
## Dublin : 8 2 : 8 1st Qu.: 13192 1st Qu.: 9.00 1:52
## EHB - Dub.: 8 3 : 8 Median : 23878 Median : 21.00
## Galway : 8 4 : 8 Mean : 34512 Mean : 36.77
## Lim. : 8 5 : 8 3rd Qu.: 46879 3rd Qu.: 52.75
## Mid HB : 8 6 : 8 Max. :210422 Max. :198.00
## (Other) :56 (Other):56
## age smr expected
## 1:26 Min. : 27.20 Min. : 1.529
## 2:26 1st Qu.: 90.08 1st Qu.: 9.296
## 3:26 Median :109.40 Median : 17.983
## 4:26 Mean :112.88 Mean : 33.647
## 3rd Qu.:134.60 3rd Qu.: 45.291
## Max. :252.30 Max. :238.267
## Complete Pooling Model
cpooling <- zelig(death ~ Region + sex, data = irlsuicide, model = "normal")## How to cite this model in Zelig:
## Kosuke Imai, Gary King, and Oliva Lau. 2007. "normal: Normal Regression for Continuous Dependent Variables" in Kosuke Imai, Gary King, and Olivia Lau, "Zelig: Everyone's Statistical Software," http://gking.harvard.edu/zeligstargazer(cpooling, type = "text", single.row = T)##
## =============================================
## Dependent variable:
## ---------------------------
## death
## ---------------------------------------------
## RegionDublin 37.000*** (12.339)
## RegionEHB - Dub. 51.750*** (12.339)
## RegionGalway -17.250 (12.339)
## RegionLim. -15.250 (12.339)
## RegionMid HB 3.000 (12.339)
## RegionMWHB - Lim. 14.500 (12.339)
## RegionNEHB 17.750 (12.339)
## RegionNWHB 4.625 (12.339)
## RegionSEHB - Wat. 34.125*** (12.339)
## RegionSHB - Cork 40.125*** (12.339)
## RegionWaterf. -19.125 (12.339)
## RegionWHB - Gal. 19.625 (12.339)
## sex1 45.000*** (4.840)
## Constant 1.125 (9.055)
## ---------------------------------------------
## Observations 104
## Log Likelihood -474.470
## Akaike Inf. Crit. 976.939
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01cpooling1 <- zelig(death ~ Region + sex + Region:sex, data = irlsuicide, model = "normal")## How to cite this model in Zelig:
## Kosuke Imai, Gary King, and Oliva Lau. 2007. "normal: Normal Regression for Continuous Dependent Variables" in Kosuke Imai, Gary King, and Olivia Lau, "Zelig: Everyone's Statistical Software," http://gking.harvard.edu/zeligstargazer(cpooling1, type = "text", single.row = T)##
## ==================================================
## Dependent variable:
## ---------------------------
## death
## --------------------------------------------------
## RegionDublin 20.500 (15.280)
## RegionEHB - Dub. 20.500 (15.280)
## RegionGalway -8.750 (15.280)
## RegionLim. -8.500 (15.280)
## RegionMid HB -2.250 (15.280)
## RegionMWHB - Lim. 2.250 (15.280)
## RegionNEHB 2.500 (15.280)
## RegionNWHB -2.000 (15.280)
## RegionSEHB - Wat. 8.000 (15.280)
## RegionSHB - Cork 13.000 (15.280)
## RegionWaterf. -8.750 (15.280)
## RegionWHB - Gal. 2.750 (15.280)
## sex1 24.750 (15.280)
## RegionDublin :sex1 33.000 (21.610)
## RegionEHB - Dub.:sex1 62.500*** (21.610)
## RegionGalway :sex1 -17.000 (21.610)
## RegionLim. :sex1 -13.500 (21.610)
## RegionMid HB:sex1 10.500 (21.610)
## RegionMWHB - Lim.:sex1 24.500 (21.610)
## RegionNEHB:sex1 30.500 (21.610)
## RegionNWHB:sex1 13.250 (21.610)
## RegionSEHB - Wat.:sex1 52.250** (21.610)
## RegionSHB - Cork:sex1 54.250** (21.610)
## RegionWaterf.:sex1 -20.750 (21.610)
## RegionWHB - Gal.:sex1 33.750 (21.610)
## Constant 11.250 (10.805)
## --------------------------------------------------
## Observations 104
## Log Likelihood -453.217
## Akaike Inf. Crit. 958.434
## ==================================================
## Note: *p<0.1; **p<0.05; ***p<0.01The first complete pooling model tells us that there is a higher rate of female suicides in EHB-Dublin (p=.01), SEHB- Wat (p=.01), SHB-Cork (p=.01), and Dublin (p=.01). The second complete pooling model takes into account the interaction between gender and geographic region. The results tell us that there is a higher rate of suicide in women in EHB-Dublin (p=.01), SEHB-Wat (p=.05), and SHB-Cork (p=.05).
No-pooling Model
dcoef <- irlsuicide %>%
group_by(Region) %>%
do(mod = lm(death ~ sex, data = .))
coef <- dcoef %>% do(data.frame(sexc = coef(.$mod)[2]))
hist(coef$sexc)
coef <- data.frame(coef)
Desc(coef)##
## -------------------------------------------------------------------------
## 'data.frame': 13 obs. of 1 variable:
## 1 $ sexc: num 24.75 57.75 87.25 7.75 11.25 ...
##
## -------------------------------------------------------------------------
## 1 - sexc (numeric)
##
## length n NAs unique 0s mean meanSE
## 13 13 0 = n 0 45 7.677
##
## .05 .10 .25 median .75 .90 .95
## 6.250 8.450 24.750 49.250 58.500 78.600 82.300
##
## rng sd vcoef mad IQR skew kurt
## 83.250 27.681 0.615 36.324 33.750 -0.056 -1.445
##
## lowest : 4, 7.75, 11.25, 24.75, 35.25
## highest: 57.75, 58.5, 77, 79, 87.25
##
## Shapiro-Wilks normality test p.value : 0.55554*Multi-level Models
m1_lme <- lmer(death ~ sex + (1|Region), data = irlsuicide)
summary(m1_lme)## Linear mixed model fit by REML ['lmerMod']
## Formula: death ~ sex + (1 | Region)
## Data: irlsuicide
##
## REML criterion at convergence: 974.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3654 -0.6604 -0.0931 0.5213 4.2800
##
## Random effects:
## Groups Name Variance Std.Dev.
## Region (Intercept) 458.2 21.41
## Residual 609.0 24.68
## Number of obs: 104, groups: Region, 13
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.269 6.853 2.082
## sex1 45.000 4.840 9.298
##
## Correlation of Fixed Effects:
## (Intr)
## sex1 -0.353m2_lme <- lmer(death ~ sex + (sex|Region), data = irlsuicide)
summary(m2_lme)## Linear mixed model fit by REML ['lmerMod']
## Formula: death ~ sex + (sex | Region)
## Data: irlsuicide
##
## REML criterion at convergence: 943.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0032 -0.2889 -0.0697 0.2623 4.0805
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Region (Intercept) 82.27 9.07
## sex1 693.73 26.34 1.00
## Residual 411.59 20.29
## Number of obs: 104, groups: Region, 13
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.269 3.774 3.781
## sex1 45.000 8.318 5.410
##
## Correlation of Fixed Effects:
## (Intr)
## sex1 0.333stargazer(cpooling, cpooling1, m1_lme, m2_lme, type = "text")##
## ==============================================================
## Dependent variable:
## ---------------------------------------
## death
## normal linear
## mixed-effects
## (1) (2) (3) (4)
## --------------------------------------------------------------
## RegionDublin 37.000*** 20.500
## (12.339) (15.280)
##
## RegionEHB - Dub. 51.750*** 20.500
## (12.339) (15.280)
##
## RegionGalway -17.250 -8.750
## (12.339) (15.280)
##
## RegionLim. -15.250 -8.500
## (12.339) (15.280)
##
## RegionMid HB 3.000 -2.250
## (12.339) (15.280)
##
## RegionMWHB - Lim. 14.500 2.250
## (12.339) (15.280)
##
## RegionNEHB 17.750 2.500
## (12.339) (15.280)
##
## RegionNWHB 4.625 -2.000
## (12.339) (15.280)
##
## RegionSEHB - Wat. 34.125*** 8.000
## (12.339) (15.280)
##
## RegionSHB - Cork 40.125*** 13.000
## (12.339) (15.280)
##
## RegionWaterf. -19.125 -8.750
## (12.339) (15.280)
##
## RegionWHB - Gal. 19.625 2.750
## (12.339) (15.280)
##
## sex1 45.000*** 24.750 45.000*** 45.000***
## (4.840) (15.280) (4.840) (8.318)
##
## RegionDublin :sex1 33.000
## (21.610)
##
## RegionEHB - Dub.:sex1 62.500***
## (21.610)
##
## RegionGalway :sex1 -17.000
## (21.610)
##
## RegionLim. :sex1 -13.500
## (21.610)
##
## RegionMid HB:sex1 10.500
## (21.610)
##
## RegionMWHB - Lim.:sex1 24.500
## (21.610)
##
## RegionNEHB:sex1 30.500
## (21.610)
##
## RegionNWHB:sex1 13.250
## (21.610)
##
## RegionSEHB - Wat.:sex1 52.250**
## (21.610)
##
## RegionSHB - Cork:sex1 54.250**
## (21.610)
##
## RegionWaterf.:sex1 -20.750
## (21.610)
##
## RegionWHB - Gal.:sex1 33.750
## (21.610)
##
## Constant 1.125 11.250 14.269** 14.269***
## (9.055) (10.805) (6.853) (3.774)
##
## --------------------------------------------------------------
## Observations 104 104 104 104
## Log Likelihood -474.470 -453.217 -487.381 -471.532
## Akaike Inf. Crit. 976.939 958.434 982.762 955.063
## Bayesian Inf. Crit. 993.340 970.930
## ==============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01library(texreg)
screenreg(list(m1_lme, m2_lme))##
## ======================================================
## Model 1 Model 2
## ------------------------------------------------------
## (Intercept) 14.27 * 14.27 ***
## (6.85) (3.77)
## sex1 45.00 *** 45.00 ***
## (4.84) (8.32)
## ------------------------------------------------------
## AIC 982.76 955.06
## BIC 993.34 970.93
## Log Likelihood -487.38 -471.53
## Num. obs. 104 104
## Num. groups: Region 13 13
## Variance: Region.(Intercept) 458.23 82.27
## Variance: Residual 609.04 411.59
## Variance: Region.sex1 693.73
## ======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05m3_lme <- lmer(death ~ sex + (sex|Region) + pop, data = irlsuicide)## Warning: Some predictor variables are on very different scales: consider
## rescaling## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.141106
## (tol = 0.002, component 2)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvaluesm4_lme <- lmer(death ~ sex*pop + (sex|Region), data = irlsuicide)## Warning: Some predictor variables are on very different scales: consider
## rescalingscreenreg(list(m3_lme, m4_lme))##
## ======================================================
## Model 1 Model 2
## ------------------------------------------------------
## (Intercept) -1.51 9.67 ***
## (2.82) (2.84)
## sex1 45.14 *** 20.69 ***
## (7.60) (4.77)
## pop 0.00 *** 0.00 **
## (0.00) (0.00)
## sex1:pop 0.00 ***
## (0.00)
## ------------------------------------------------------
## AIC 922.69 883.69
## BIC 941.20 904.85
## Log Likelihood -454.34 -433.85
## Num. obs. 104 104
## Num. groups: Region 13 13
## Variance: Region.(Intercept) 0.00 29.06
## Variance: Region.sex1 623.41 144.48
## Variance: Residual 256.70 147.01
## ======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05Centering the Independent Variables
library(magrittr)
irlsuicide %<>% mutate(cpop = pop - mean(pop))m4a_lme <- lmer(death ~ sex*cpop + (sex|Region), data = irlsuicide)## Warning: Some predictor variables are on very different scales: consider
## rescalingscreenreg(list(m4_lme, m4a_lme))##
## ======================================================
## Model 1 Model 2
## ------------------------------------------------------
## (Intercept) 9.67 *** 14.25 ***
## (2.84) (2.25)
## sex1 20.69 *** 45.15 ***
## (4.77) (4.09)
## pop 0.00 **
## (0.00)
## sex1:pop 0.00 ***
## (0.00)
## cpop 0.00 **
## (0.00)
## sex1:cpop 0.00 ***
## (0.00)
## ------------------------------------------------------
## AIC 883.69 883.69
## BIC 904.85 904.85
## Log Likelihood -433.85 -433.85
## Num. obs. 104 104
## Num. groups: Region 13 13
## Variance: Region.(Intercept) 29.06 29.06
## Variance: Region.sex1 144.48 144.48
## Variance: Residual 147.01 147.01
## ======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05Intra-Class Correlation
m0_lme <- lmer(death ~ 1 + (1|Region), data = irlsuicide)
summary(m0_lme)## Linear mixed model fit by REML ['lmerMod']
## Formula: death ~ 1 + (1 | Region)
## Data: irlsuicide
##
## REML criterion at convergence: 1041
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3181 -0.6837 -0.3082 0.3124 3.8787
##
## Random effects:
## Groups Name Variance Std.Dev.
## Region (Intercept) 386.7 19.67
## Residual 1180.9 34.36
## Number of obs: 104, groups: Region, 13
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 36.769 6.411 5.735To be honest, I wasn’t sure where to find the information to help me interpret these results. I decided to continue with the Bristol data as well, in case that is easier to interpret.