setwd("/media/guerzhoy/Windows/STA303/lectures/lec10")
#setwd("c:/STA303/lectures/lec10")
library(lme4)
library(foreign)
Let’s fit the model:
\[P(y_i) = logistic(\alpha^{state}_{j[i]}+\beta_{female}I_{female, i}+\beta_{black}I_{black, i})\] \[\alpha^{state}_j\sim N(\mu_{\alpha}, \sigma^{2}_{state})\]
polls <- read.dta ("polls.dta") #The 1988 Election: Bush (Rep.) vs Dukakis (Dem.)
fit <- glmer (bush ~ black + female + (1 | state), family=binomial(link="logit"), data=polls)
summary(fit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: bush ~ black + female + (1 | state)
## Data: polls
##
## AIC BIC logLik deviance df.resid
## 15153.2 15182.6 -7572.6 15145.2 11562
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8924 -1.0901 0.6273 0.8334 2.6628
##
## Random effects:
## Groups Name Variance Std.Dev.
## state (Intercept) 0.1677 0.4096
## Number of obs: 11566, groups: state, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.43318 0.06946 6.236 4.48e-10 ***
## black -1.81573 0.08759 -20.729 < 2e-16 ***
## female -0.11557 0.03936 -2.936 0.00332 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) black
## black -0.058
## female -0.325 -0.026
ranef(fit)
## $state
## (Intercept)
## 1 0.701039244
## 3 0.074583311
## 4 0.016954616
## 5 0.004885914
## 6 0.068798464
## 7 -0.163072325
## 8 -0.400761663
## 9 -0.548664395
## 10 0.290717734
## 11 0.207809454
## 13 -0.216403739
## 14 -0.061265498
## 15 0.402556668
## 16 -0.649097333
## 17 0.407642723
## 18 0.063350658
## 19 0.639696387
## 20 -0.222406251
## 21 -0.090855219
## 22 -0.460632157
## 23 0.045245224
## 24 -0.260655157
## 25 0.842507237
## 26 -0.301495690
## 27 -0.460279300
## 28 -0.015754855
## 29 0.070918191
## 30 0.142933031
## 31 -0.048600192
## 32 -0.142062118
## 33 -0.456429512
## 34 0.243043099
## 35 -0.068780655
## 36 0.180568983
## 37 -0.067709578
## 38 -0.248834354
## 39 -0.186890074
## 40 -0.782204047
## 41 0.617868187
## 42 -0.050871617
## 43 0.534786714
## 44 0.059228653
## 45 0.612534868
## 46 -0.009517583
## 47 0.615168926
## 48 -0.427131952
## 49 -0.273427716
## 50 -0.305426771
## 51 0.010000341
It is better to include an interaction so that we can estimate the probability of voting Republican for each demographic, broken down by sex, race, and age (age is categorical as well).
We are now fitting the model
\[P(y_i) = logistic(\alpha^{state}_{j[i]}+\alpha^{age}_{k[i]}+\beta_{female}I_{female, i}+\beta_{black}I_{black, i})+\beta{female, black}I_{female, i}I_{black, i} \] \[\alpha_j\sim N(\mu^{state}_{\alpha}, \sigma^{2}_{state})\] \[\alpha_k\sim N(\mu^{age}_{\alpha}, \sigma^{2}_{state})\]
polls$age <- factor(polls$age)
fit2 <- glmer (bush ~ black + female + black:female + (1 | state) + (1|age), family=binomial(link="logit"), data=polls)
summary(fit2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: bush ~ black + female + black:female + (1 | state) + (1 | age)
## Data: polls
##
## AIC BIC logLik deviance df.resid
## 15140.7 15184.9 -7564.4 15128.7 11560
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9114 -1.0874 0.6369 0.8362 2.9324
##
## Random effects:
## Groups Name Variance Std.Dev.
## state (Intercept) 0.170235 0.41260
## age (Intercept) 0.009617 0.09807
## Number of obs: 11566, groups: state, 49; age, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44214 0.08577 5.155 2.54e-07 ***
## black -2.07706 0.14985 -13.861 < 2e-16 ***
## female -0.12832 0.04058 -3.162 0.00156 **
## black:female 0.36579 0.18143 2.016 0.04379 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) black female
## black -0.073
## female -0.274 0.154
## black:femal 0.061 -0.809 -0.225
ranef(fit2)
## $state
## (Intercept)
## 1 0.697598365
## 3 0.088795779
## 4 0.033468238
## 5 -0.005117006
## 6 0.058133310
## 7 -0.171028036
## 8 -0.393215992
## 9 -0.562590848
## 10 0.297023318
## 11 0.201312360
## 13 -0.221526511
## 14 -0.056518189
## 15 0.399951907
## 16 -0.659906543
## 17 0.408706381
## 18 0.068361644
## 19 0.662854721
## 20 -0.227841658
## 21 -0.094300270
## 22 -0.468486535
## 23 0.043231963
## 24 -0.265370155
## 25 0.843530393
## 26 -0.298770127
## 27 -0.456988221
## 28 -0.012490469
## 29 0.095505637
## 30 0.148479686
## 31 -0.058204772
## 32 -0.137222851
## 33 -0.457232223
## 34 0.255207475
## 35 -0.072442629
## 36 0.179632905
## 37 -0.060814394
## 38 -0.249012926
## 39 -0.179461325
## 40 -0.783565322
## 41 0.625708770
## 42 -0.054728945
## 43 0.541487692
## 44 0.058125710
## 45 0.600242505
## 46 -0.010054744
## 47 0.597783221
## 48 -0.422880005
## 49 -0.269559958
## 50 -0.315014257
## 51 0.014991172
##
## $age
## (Intercept)
## 1 0.138306962
## 2 -0.009160998
## 3 -0.059495324
## 4 -0.072148251