Homework 9
Robert Johnson
Sunday, April 26, 2015
## 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.
## ####
## 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=stargazer
##
##
## 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, union
##
## Loading required package: Matrix
## Loading required package: Rcpp
##
## Attaching package: 'npmlreg'
##
## The following object is masked from 'package:Matrix':
##
## expand## 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/zelig##
## =============================================
## Dependent variable:
## ---------------------------
## death
## ---------------------------------------------
## sex1 45.000*** (6.326)
## Constant 14.269*** (4.473)
## ---------------------------------------------
## Observations 104
## Log Likelihood -508.819
## Akaike Inf. Crit. 1,021.637
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
##
## -------------------------------------------------------------------------
## '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.55554The data I am using is a study on Suicide Rates in the Republic of Ireland 1989-1998.
Region a factor with levels Cork , Dublin , EHB - Dub., Galway, Lim., Mid HB, MWHB-Lim., NEHB, NWHB, SEHB-Wat., SHB-Cork, Waterf., WHB-Gal.. sex a factor for gender with levels 0 (female) and 1 (male). death a numeric vector giving the total number of deaths. expected a numeric vector with `expected’ number of cases obtained from a reference population (Ahlbom, 1993). (death in suicides I assume)
## 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.353## 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.333##
## =================================================
## Dependent variable:
## -----------------------------
## death
## normal linear
## mixed-effects
## (1) (2) (3)
## -------------------------------------------------
## sex1 45.000*** 45.000*** 45.000***
## (6.326) (4.840) (8.318)
##
## Constant 14.269*** 14.269** 14.269***
## (4.473) (6.853) (3.774)
##
## -------------------------------------------------
## Observations 104 104 104
## Log Likelihood -508.819 -487.381 -471.532
## Akaike Inf. Crit. 1,021.637 982.762 955.063
## Bayesian Inf. Crit. 993.340 970.930
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01## Version: 1.34
## Date: 2014-10-31
## Author: Philip Leifeld (University of Konstanz)
##
## Please cite the JSS article in your publications -- see citation("texreg").In the above model we see that all three models are statistically significant at a 99% confidence level. To see which model to use we look at the Akaike, model 1 has the highest rate followed by 2 then 3.
The above shows that there is a statistical correlation between sex and death by suicide, specifically sex1 (male). With a p value of .01 or lower meaning that it is at a 99% confidence level.
The pooling model puts a median of 49.250 because of the large varience in range of deaths by sucidie.
The table below shows that there is a significant between death, sex and region.
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) + expected, data = irlsuicide)
m4_lme <- lmer(age ~ sex*expected + (sex|Region), data = irlsuicide)
screenreg(list(m3_lme, m4_lme))##
## ======================================================
## Model 1 Model 2
## ------------------------------------------------------
## (Intercept) 1.05 2.45 ***
## (1.35) (0.24)
## sex1 5.44 0.47
## (3.83) (0.34)
## expected 0.98 *** 0.00
## (0.04) (0.01)
## sex1:expected -0.01
## (0.01)
## ------------------------------------------------------
## AIC 763.59 349.03
## BIC 782.10 370.19
## Log Likelihood -374.79 -166.52
## Num. obs. 104 104
## Num. groups: Region 13 13
## Variance: Region.(Intercept) 5.05 0.00
## Variance: Region.sex1 127.60 0.00
## Variance: Residual 59.80 1.23
## ======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05library(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:texreg':
##
## extractirlsuicide %<>% mutate(death = expected - mean(expected))
m4a_lme <- lmer(death ~ sex*expected + (sex|Region), data = irlsuicide)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?screenreg(list(m4_lme, m4a_lme))##
## =======================================================
## Model 1 Model 2
## -------------------------------------------------------
## (Intercept) 2.45 *** -33.65 ***
## (0.24) (0.00)
## sex1 0.47 0.00 ***
## (0.34) (0.00)
## expected 0.00 1.00 ***
## (0.01) (0.00)
## sex1:expected -0.01 -0.00 ***
## (0.01) (0.00)
## -------------------------------------------------------
## AIC 349.03 -5808.35
## BIC 370.19 -5787.20
## Log Likelihood -166.52 2912.18
## Num. obs. 104 104
## Num. groups: Region 13 13
## Variance: Region.(Intercept) 0.00 0.00
## Variance: Region.sex1 0.00 0.00
## Variance: Residual 1.23 0.00
## =======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05m0_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: 1032
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8822 -0.5660 -0.2454 0.2597 4.8004
##
## Random effects:
## Groups Name Variance Std.Dev.
## Region (Intercept) 506.4 22.50
## Residual 1045.4 32.33
## Number of obs: 104, groups: Region, 13
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -7.494e-15 7.000e+00 0