DEM 7283 - Multi-level Models Assignment
Paulina Cano
April 2, 2018
All the variables in this dataset were compiled the Behavioral Risk Factor Surveillance System (BRFFS).
For this assignment, we include as a continous variable the average alcoholic drinks consumed per day in past month (30 days).For individual predictors, we look at the level of education and race/ethnicity for all female participants.
Our higher level variable is the Metropolitan Metro areas included in the same dataset.
#load brfss
library(car)
library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.1. https://CRAN.R-project.org/package=stargazerlibrary(survey)## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartlibrary(sjPlot)
library(ggplot2)
library(pander)
library(knitr)load("C:/Users/PCMcC/Documents/Statistics II/Data for project/brfss16_mmsa.Rdata")
set.seed(12345)
#samps<-sample(1:nrow(brfss16m), size = 249011, replace=F)
#brfss16m<-brfss16m[samps,]
#The names in the data are very ugly, so I make them less ugly
nams<-names(brfss16m)
#we see some names are lower case, some are upper and some have a little _ in the first position. This is a nightmare.
newnames<-gsub(pattern = "_",replacement = "",x = nams)
names(brfss16m)<-tolower(newnames)
nams<-names(brfss16m)
#we see some names are lower case, some are upper and some have a little _ in the first position. This is a nightmare.
newnames<-gsub(pattern = "_",replacement = "",x = nams)
names(brfss16m)<-tolower(newnames)Recode of Continuous Variable
#Avg alcoholic drinks per day in past 30.
#Survey Question: One drink is equivalent to a 12-ounce beer, a 5-ounce glass of wine, or a drink with one shot of liquor. During the past 30 days, on the days when you drank, about how many drinks did you drink on the average? (A 40 ounce beer would count as 3 drinks, or a cocktail drink with 2 shots would count as 2 drinks.)
brfss16m$daysdrink<-recode(brfss16m$avedrnk2, recodes = "77=NA; 99=NA")
summary(brfss16m$daysdrink)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 1.00 2.00 2.17 2.00 98.00 122698Recode of other variables
#sex
brfss16m$female<-ifelse(brfss16m$sex==2, 2, 0)
#BMI
brfss16m$bmi<-ifelse(is.na(brfss16m$bmi5)==T, NA, brfss16m$bmi5/100)
#Healthy days
brfss16m$healthdays<-recode(brfss16m$physhlth, recodes = "88=0; 77=NA; 99=NA")
#Healthy mental health days
brfss16m$healthmdays<-recode(brfss16m$menthlth, recodes = "88=0; 77=NA; 99=NA")
#Education
brfss16m$education<-recode(brfss16m$educa, recodes="1='no school'; 2='Grade 1-8'; 3='Grade 9-11';4='high school'; 5='some college'; 6='college or more'; 9='NA'",as.factor.result = T)
brfss16m$race_eth<-relevel(brfss16m$education, ref = "college or more")
brfss16m$badhealth<-recode(brfss16m$genhlth, recodes="4:5=1; 1:3=0; else=NA")
#race/ethnicity
brfss16m$black<-recode(brfss16m$racegr3, recodes="2=1; 9=NA; else=0")
brfss16m$white<-recode(brfss16m$racegr3, recodes="1=1; 9=NA; else=0")
brfss16m$other<-recode(brfss16m$racegr3, recodes="3:4=1; 9=NA; else=0")
brfss16m$hispanic<-recode(brfss16m$racegr3, recodes="5=1; 9=NA; else=0")
brfss16m$race_eth<-recode(brfss16m$racegr3,
recodes="1='nhwhite'; 2='nh black'; 3='nh other';4='nh multirace'; 5='hispanic'; else=NA",
as.factor.result = T)
brfss16m$race_eth<-relevel(brfss16m$race_eth, ref = "nhwhite")
#insurance
brfss16m$ins<-recode(brfss16m$hlthpln1, recodes ="7:9=NA; 1=1;2=0")
#income grouping
brfss16m$inc<-recode(brfss16m$incomg, recodes = "9= NA;1='1_lt15k'; 2='2_15-25k';3='3_25-35k';4='4_35-50k';5='5_50kplus'", as.factor.result = T)
brfss16m$inc<-as.ordered(brfss16m$inc)
#education level
brfss16m$educ<-recode(brfss16m$educa,
recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA",
as.factor.result=T)
brfss16m$educ<-relevel(brfss16m$educ, ref='2hsgrad')
#employment
brfss16m$employ<-recode(brfss16m$employ1,
recodes="1:2='employloyed'; 2:6='nilf'; 7='retired'; 8='unable'; else=NA",
as.factor.result=T)
brfss16m$employ<-relevel(brfss16m$employ, ref='employloyed')
#marital status
brfss16m$marst<-recode(brfss16m$marital,
recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA",
as.factor.result=T)
brfss16m$marst<-relevel(brfss16m$marst, ref='married')
#Age cut into intervals
brfss16m$agec<-cut(brfss16m$age80, breaks=c(0,24,39,59,79,99))
#BMI, in the brfss16ma the bmi variable has 2 implied decimal places,
#so we must divide by 100 to get real bmi's
brfss16m$bmi<-brfss16m$bmi5/100
#smoking currently
brfss16m$smoke<-recode(brfss16m$smoker3,
recodes="1:2=1; 3:4=0; else=NA")
#brfss16m$smoke<-relevel(brfss16m$smoke, ref = "NeverSmoked")library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnewdata<-brfss16m%>%
select(education,race_eth,daysdrink,sex,female,mmsaname, mmsa)%>%
filter(complete.cases(.))Information on the MSAs
There information on MSA shows that there are a total of 143 metropolitan areas.
head(data.frame(name=table(brfss16m$mmsaname),id=unique(brfss16m$mmsa)))## name.Var1
## 1 Albany-Schenectady-Troy, NY, Metropolitan Statistical Area
## 2 Albuquerque, NM, Metropolitan Statistical Area
## 3 Allentown-Bethlehem-Easton, PA-NJ, Metropolitan Statistical Area
## 4 Anchorage, AK, Metropolitan Statistical Area
## 5 Atlanta-Sandy Springs-Roswell, GA, Metropolitan Statistical Area
## 6 Augusta-Richmond County, GA-SC, Metropolitan Statistical Area
## name.Freq id
## 1 2901 10580
## 2 1300 10740
## 3 675 10900
## 4 1028 11260
## 5 2444 12060
## 6 849 12260length(table(brfss16m$mmsa))## [1] 143Basic hierarchical model for means of each MSA:
library(lme4)
library(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(arm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## arm (Version 1.9-3, built: 2016-11-21)## Working directory is C:/Users/PCMcC/Documents/Statistics II/Data for project/highered##
## Attaching package: 'arm'## The following object is masked from 'package:car':
##
## logit##Scaling continous variable
newdata$daysdrinkz<-as.numeric(scale(newdata$daysdrink, center=T, scale=T))
fit<-lmer(daysdrinkz~(1|mmsa),data=newdata)
arm::display(fit, detail=T)## lme4::lmer(formula = daysdrinkz ~ (1 | mmsa), data = newdata)
## coef.est coef.se t value
## (Intercept) 0.02 0.01 1.43
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.15
## Residual 0.99
## ---
## number of obs: 124408, groups: mmsa, 143
## AIC = 350919, DIC = 350899.5
## deviance = 350906.4head(ranef(fit))## $mmsa
## (Intercept)
## 10580 -0.008568762
## 10740 -0.077010741
## 10900 0.035958396
## 11260 -0.059864846
## 12060 -0.007580920
## 12260 -0.019246920
## 12420 -0.080437432
## 12580 -0.109983518
## 12940 -0.039243533
## 13220 0.001236381
## 13620 -0.148272105
## 13740 -0.008435657
## 13780 0.028505947
## 13820 0.059263503
## 13900 0.014638366
## 14260 -0.046403481
## 14454 0.054502622
## 15380 -0.017678933
## 15540 -0.124566952
## 15764 -0.048072657
## 15804 -0.065334804
## 16300 -0.026626011
## 16620 0.098393661
## 16700 -0.013314461
## 16740 0.043665756
## 16860 -0.081342785
## 16980 0.036404862
## 17140 0.151003764
## 17200 -0.154574336
## 17460 -0.019388038
## 17780 -0.134136059
## 17820 -0.134080688
## 17900 -0.004047075
## 18140 -0.020694479
## 18580 -0.007781237
## 18880 0.012685502
## 19060 0.019592055
## 19124 -0.030128550
## 19380 -0.080523968
## 19660 0.007179178
## 19740 -0.102981269
## 19780 -0.030492996
## 20260 -0.007271288
## 20524 -0.088121644
## 21340 0.191116085
## 22020 -0.021418547
## 22220 -0.065418798
## 23060 -0.047954368
## 23104 -0.063116580
## 23540 0.030856565
## 24020 0.036408452
## 24220 0.026072636
## 24260 -0.045736167
## 24340 0.067812757
## 24860 -0.046428904
## 25180 -0.001850287
## 25540 -0.152043938
## 25940 -0.149800361
## 26420 -0.030034942
## 26580 0.124344982
## 26900 -0.072308353
## 27140 0.053349473
## 27260 0.050371542
## 28140 -0.057202822
## 28700 -0.016893799
## 28940 0.061152978
## 29620 0.059697404
## 30700 0.023496453
## 30780 -0.088205716
## 30860 0.034345263
## 31080 -0.008373256
## 31140 0.068489304
## 32820 -0.040119322
## 33100 0.108206180
## 33340 -0.025843251
## 33460 -0.049655543
## 33500 -0.014039228
## 33874 -0.059994134
## 34820 -0.022592960
## 34980 0.008417130
## 35004 -0.079674269
## 35084 -0.083708897
## 35380 0.031842769
## 35614 0.004428719
## 35740 0.007857725
## 35820 -0.020310638
## 35840 -0.058813018
## 36084 -0.022438685
## 36260 0.140847781
## 36420 -0.056410076
## 36540 0.052069754
## 36740 0.098639490
## 37460 0.131400410
## 37860 0.038669240
## 37964 0.070041254
## 38060 -0.118070227
## 38300 0.068884346
## 38860 -0.142466567
## 38900 -0.067145923
## 38940 -0.015841932
## 39300 -0.084315566
## 39340 0.407580426
## 39580 0.002651405
## 39660 -0.149468960
## 39900 0.038437489
## 40060 -0.032507645
## 40140 0.283020161
## 40340 -0.015360563
## 40380 -0.060423914
## 40484 -0.081196993
## 40900 -0.004116998
## 41060 0.054976826
## 41180 0.034098932
## 41420 -0.054037585
## 41540 -0.020630292
## 41620 0.204287243
## 41700 0.069126573
## 41884 -0.033016403
## 41940 -0.074313440
## 41980 1.335056672
## 42420 -0.050370856
## 42644 -0.144660090
## 43524 -0.176217960
## 43580 0.016460400
## 43620 -0.107936920
## 43900 0.013758756
## 44060 -0.100048606
## 44140 -0.031301019
## 45060 0.033491726
## 45220 0.009860103
## 45300 0.093256707
## 45780 0.010798566
## 45820 0.032546210
## 46140 -0.071580239
## 46220 0.133499205
## 46540 0.030793838
## 47260 0.015543617
## 47664 -0.002503177
## 47894 -0.044396752
## 48620 -0.034659453
## 48660 -0.063385146
## 48864 0.008136317
## 49340 -0.002663623dim(ranef(fit)$mmsa)## [1] 143 1rate<- fixef(fit)+ranef(fit)$mmsa
est<-data.frame(rate =rate, id=rownames(ranef(fit)$mmsa))meandaysdrink<-mean(newdata$daysdrink, na.rm=T)
sddaysdrink<-sd(newdata$daysdrink, na.rm=T)newdata$preddaysdrink<-sddaysdrink*(fitted(fit))+meandaysdrink
citymeans<-aggregate(cbind(daysdrink,preddaysdrink)~mmsaname,newdata,mean)
head(citymeans, n=20) ## mmsaname
## 1 Albany-Schenectady-Troy, NY, Metropolitan Statistical Area
## 2 Albuquerque, NM, Metropolitan Statistical Area
## 3 Allentown-Bethlehem-Easton, PA-NJ, Metropolitan Statistical Area
## 4 Anchorage, AK, Metropolitan Statistical Area
## 5 Atlanta-Sandy Springs-Roswell, GA, Metropolitan Statistical Area
## 6 Augusta-Richmond County, GA-SC, Metropolitan Statistical Area
## 7 Austin-Round Rock, TX, Metropolitan Statistical Area
## 8 Baltimore-Columbia-Towson, MD, Metropolitan Statistical Area
## 9 Baton Rouge, LA, Metropolitan Statistical Area
## 10 Beckley, WV, Metropolitan Statistical Area
## 11 Berlin, NH-VT, Micropolitan Statistical Area
## 12 Billings, MT, Metropolitan Statistical Area
## 13 Binghamton, NY, Metropolitan Statistical Area
## 14 Birmingham-Hoover, AL, Metropolitan Statistical Area
## 15 Bismarck, ND, Metropolitan Statistical Area
## 16 Boise City, ID, Metropolitan Statistical Area
## 17 Boston, MA, Metropolitan Division
## 18 Buffalo-Cheektowaga-Niagara Falls, NY, Metropolitan Statistical Area
## 19 Burlington-South Burlington, VT, Metropolitan Statistical Area
## 20 Cambridge-Newton-Framingham, MA, Metropolitan Division
## daysdrink preddaysdrink
## 1 2.189342 2.189934
## 2 2.012678 2.026125
## 3 2.307479 2.296505
## 4 2.054717 2.067162
## 5 2.191547 2.192298
## 6 2.158854 2.164377
## 7 2.007910 2.017924
## 8 1.943785 1.947208
## 9 2.101399 2.116517
## 10 2.214286 2.213401
## 11 1.787500 1.855569
## 12 2.187302 2.190252
## 13 2.285408 2.278668
## 14 2.366890 2.352283
## 15 2.248025 2.245477
## 16 2.091200 2.099380
## 17 2.345427 2.340888
## 18 2.165123 2.168129
## 19 1.898023 1.912305
## 20 2.092253 2.095385A summary of our citymeans object shows that the max number of drinks consumed is 5.5 and the minimum 1.77 among the MSA’s.
summary(citymeans)## mmsaname daysdrink preddaysdrink
## Length:143 Min. :1.778 Min. :1.789
## Class :character 1st Qu.:2.057 1st Qu.:2.067
## Mode :character Median :2.169 Median :2.174
## Mean :2.212 Mean :2.210
## 3rd Qu.:2.295 3rd Qu.:2.289
## Max. :5.572 Max. :5.406ANOVA
After computing an analysis of variance test, we see that there is a small F test produced which shows that there is a relationship between our variables.
fit.an<-lm(daysdrinkz~as.factor(mmsa), newdata)
anova(fit.an)## Analysis of Variance Table
##
## Response: daysdrinkz
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(mmsa) 142 2644 18.6189 19.001 < 2.2e-16 ***
## Residuals 124265 121763 0.9799
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Basic Hierarchical Models
fit<-lmer (daysdrinkz~(1|mmsa), data=newdata)
arm::display(fit, detail=T)## lme4::lmer(formula = daysdrinkz ~ (1 | mmsa), data = newdata)
## coef.est coef.se t value
## (Intercept) 0.02 0.01 1.43
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.15
## Residual 0.99
## ---
## number of obs: 124408, groups: mmsa, 143
## AIC = 350919, DIC = 350899.5
## deviance = 350906.4Model with Individual Predictor:Female Gender and Education
After including the variable of education for all female survey respondents, we observe an increasing positive relationship with less levels of education when compared to females with a college education or more. In terms of race/ethnicity, we see that when compared to the white population, females who are hispanic or multirace are more likely than white females to report having consumed more alcoholic drinks than white females and black.
fit2<-lmer(daysdrinkz~female+education+race_eth+(1|mmsa), data=newdata, na.action=na.omit)
arm::display(fit2, detail=T)## lme4::lmer(formula = daysdrinkz ~ female + education + race_eth +
## (1 | mmsa), data = newdata, na.action = na.omit)
## coef.est coef.se t value
## (Intercept) 0.05 0.01 4.12
## female -0.16 0.00 -57.43
## educationGrade 1-8 0.50 0.03 17.42
## educationGrade 9-11 0.44 0.02 24.71
## educationhigh school 0.24 0.01 31.97
## educationNA 0.00 0.08 0.03
## educationno school 0.38 0.11 3.40
## educationsome college 0.15 0.01 22.37
## race_ethhispanic 0.27 0.01 22.70
## race_ethnh black 0.01 0.01 0.64
## race_ethnh multirace 0.14 0.02 6.06
## race_ethnh other 0.06 0.02 3.64
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.12
## Residual 0.97
## ---
## number of obs: 124408, groups: mmsa, 143
## AIC = 344836, DIC = 344653.3
## deviance = 344730.6In doing a likelihood ratio test to see if our individual level variables are explaining anything about our outcome. We observe a significant effect with the various levels of education and race/ethnicity.
anova(fit, fit2)## refitting model(s) with ML (instead of REML)## Data: newdata
## Models:
## object: daysdrinkz ~ (1 | mmsa)
## ..1: daysdrinkz ~ female + education + race_eth + (1 | mmsa)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 3 350912 350942 -175453 350906
## ..1 14 344759 344895 -172365 344731 6175.9 11 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Random slope models
library(stargazer)
#stargazer(fit.mix3, fit.mix2, type = "html", out = "~/Google Drive/classes/dem7283/class18/lmerout.html")
library(lme4)
library(lmerTest)
library(arm)Fitting the Basic Multilevel Model
The standard deviation is .12 and the residual standard deviation .97. Less than 1% of the variance in the average number of drinks consumed is explained by the differences in MSAs.
The fixed effects continue to show that females with less levels of education are more likely to report consuming more alcohol as well as Hispanic and multirace participants.
fit.mix<-lmer(daysdrinkz~female+education+race_eth+(1|mmsa), data=newdata)
#do a test for the random effect
rand(fit.mix)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## mmsa 1282 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1display(fit.mix, detail=T)## lme4::lmer(formula = daysdrinkz ~ female + education + race_eth +
## (1 | mmsa), data = newdata)
## coef.est coef.se t value
## (Intercept) 0.05 0.01 4.12
## female -0.16 0.00 -57.43
## educationGrade 1-8 0.50 0.03 17.42
## educationGrade 9-11 0.44 0.02 24.71
## educationhigh school 0.24 0.01 31.97
## educationNA 0.00 0.08 0.03
## educationno school 0.38 0.11 3.40
## educationsome college 0.15 0.01 22.37
## race_ethhispanic 0.27 0.01 22.70
## race_ethnh black 0.01 0.01 0.64
## race_ethnh multirace 0.14 0.02 6.06
## race_ethnh other 0.06 0.02 3.64
##
## Error terms:
## Groups Name Std.Dev.
## mmsa (Intercept) 0.12
## Residual 0.97
## ---
## number of obs: 124408, groups: mmsa, 143
## AIC = 344836, DIC = 344653.3
## deviance = 344730.6Fixed Effect Model
fixcoef<-fit.an$coefficients[1]+fit.an$coefficients[-1]
fixcoef<-(fixcoef*sddaysdrink)+meandaysdrink
plot(NULL,ylim=c(1, 5), xlim=c(0,1),
ylab="Intercept", xlab="") # get the ylim from a summary of rancoefs1
title(main="Fixed Effect Model")
for (i in 1: length(fixcoef)[1]){
#I plug in a beta, here it's the effect of age from fit.mix
abline(a=fixcoef[i], b=0, lwd=1.5, col="green")
}
#Visualization of the random intercepts by plotting them
rancoefs1<-ranef(fit.mix)$mmsa+fixef(fit.mix)[1]
rancoefs1<-(rancoefs1*sddaysdrink)+meandaysdrink
summary(rancoefs1)## (Intercept)
## Min. :1.959
## 1st Qu.:2.156
## Median :2.256
## Mean :2.278
## 3rd Qu.:2.344
## Max. :4.811plot(NULL,ylim=c(1,5), xlim=c(0,1),
ylab="Intercept", xlab="female") # get the ylim from a summary of rancoefs1
title(main="Random Intercept Models")
for (i in 1: length(rancoefs1[,1])){
#I plug in a beta, here it's the effect of age from fit.mix
abline(a=rancoefs1[i,1], b=.0, lwd=1.5, col="blue")
}
###Random Slope Model
The results from the random slope model shows results that are slightly different. The females more likely to drink higher quantities of alcohol are those in grade schools 1-8 rather than those without any education. The trends in race/ethnicity continue to be the same as those observed in previous analysis. It is possible that through this model we can observe more detailed differences in the education levels than with other models.
library(stargazer)
fit.mix3<-lmer(daysdrinkz~female+education+race_eth+(female|mmsa), newdata, REML=F)
#fit.mix3<-refit(fit.mix3)
display(fit.mix3, detail=T)## lme4::lmer(formula = daysdrinkz ~ female + education + race_eth +
## (female | mmsa), data = newdata, REML = F)
## coef.est coef.se t value
## (Intercept) 0.05 0.01 3.47
## female -0.16 0.00 -33.24
## educationGrade 1-8 0.50 0.03 17.26
## educationGrade 9-11 0.44 0.02 24.53
## educationhigh school 0.23 0.01 31.76
## educationNA 0.00 0.08 0.04
## educationno school 0.38 0.11 3.36
## educationsome college 0.15 0.01 22.25
## race_ethhispanic 0.27 0.01 22.91
## race_ethnh black 0.01 0.01 0.88
## race_ethnh multirace 0.14 0.02 6.07
## race_ethnh other 0.06 0.02 3.72
##
## Error terms:
## Groups Name Std.Dev. Corr
## mmsa (Intercept) 0.16
## female 0.05 -1.00
## Residual 0.96
## ---
## number of obs: 124408, groups: mmsa, 143
## AIC = 344496, DIC = 344464.1
## deviance = 344464.1Model Comparison with a LRT
A model comparison shows that the random slope model fits better.
anova(fit.mix, fit.mix3)## refitting model(s) with ML (instead of REML)## Data: newdata
## Models:
## object: daysdrinkz ~ female + education + race_eth + (1 | mmsa)
## ..1: daysdrinkz ~ female + education + race_eth + (female | mmsa)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 14 344759 344895 -172365 344731
## ..1 16 344496 344652 -172232 344464 266.44 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#the random slope model fits better
#plot random slopes and intercepts
rancoefs2<-ranef(fit.mix3)
plot(NULL, ylim=c(-.3, 2), xlim=c(0,1),ylab="Intercept", xlab="Female Gender == 1")
title (main="Random Slope and Intercept Model - Female Gender by Metro")
cols<-rainbow(length(unique(newdata$mmsa)))
for (i in 1: dim(rancoefs2$mmsa)[1]){
abline(a=fixef(fit.mix3)["(Intercept)"]+rancoefs2$mmsa[[1]][i],
b=fixef(fit.mix3)["female"]+rancoefs2$mmsa[, "female"][i], col=cols[i])
}