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