Lecture notes

Fundamental problem

Names - Mixed effects models - Hierarchical linear models - Random effects models

These all refer to the same suite of models, but different disciplines/authors/applications have managed to confuse us all.

Perspectives in Multi-level Modeling

Basic multilevel perspectives * The Social epidemiology perspective * General ecological models * Longitudinal models for change

Each of these situate humans within some higher, contextual level that is assumed to either influence or be influenced by their actions/behaviors

And we would ideally like to be able to incorporate covariates at the individual and higher level that could influence behaviors

Multi-stage sampling

Figures from Kawachi and Berkman (2003)

Here’s a picture of this:

Multistage Sampling

Multistage Sampling

Multi-level propositions

When we have a research statement that involves individuals within some context, this is a multi-level proposition. In this sense, we are interested in questions that relate variables at different levels, the micro and the macro. This also holds in general if a sample was collected with a multi-stage sampling scheme.

In a multilevel proposition, variables are present at two different levels, and we are interested in the relationship between both the micro and macro level association with our outcome, y.

Multi-level

Multi-level

This can be contrasted with a purely micro level proposition, where all our observed variables are the level of the individual

micro-level

micro-level

Likewise, if we are only interested in the relationship between macro level variables, we have this situation:

macro-level

macro-level

Macro - micro propositions

We commonly encounter the situation where a macro level variable effects a micro level outcome. This can happen in several different ways.

macro-micro

macro-micro

The first case is a macro to micro proposition, which may be exemplified by a statement such as: “Individuals in areas with high environmental contamination leads to higher risk of death”.

Whereas the second frame illustrates a more specific special case, where there is a macro level effect, net of the individual level predictor, and may be stated “For individuals with the a given level of education, living in areas with high environmental contamination leads to higher risk of death”.

The last panel illustrates what is known as a cross level interaction, or a macro-micro interaction. This is where the relationship between x and y is dependent on Z. This leads to the statement “Individuals with low levels of education, living in areas with high environmental contamination have higher risk of death”.

Longitudinal Models

These kinds of models are used when you have observations on the same individuals over multiple time points * Don’t have to be the same times/ages for each observation + More flexibility * These models treat the individual as the higher level unit, and you are interested in studying * Change over time within an individual * Impacts of prior circumstances on later outcomes * Two modeling strategies will allow us to consider individual change over time, versus population averaged change over time.

Multilevel data Preface

Linear Mixed Model

\[y_{ij} = \mu + u_{j} + e_{ij}\]

where \[\mu\] is the grand mean, \(u_{j}\) is the fixed effect of the \(j^{th}\) group and \(e_{ij}\) is the residual for each individual

======================================================== This model assumes that you are capturing all variation in y by the group factor differences, \(u_{j}\) alone.

If you have all your groups,

and your only predictor is the group (factor) level,

and if you expect there to be directional differences across groups a priori,

then this is probably the model for you.

You might use this framework if you want to crudely model the effect of region of residence in an analysis. + i.e. is the mean different across my region?

ANOVA and ANCOVA

\[ y_{ij} = \mu + u_{j} + e_{ij} \]

\[ \mu + u_{j} \]

\[ y_{ij} = \mu + u_{j}*\beta x_{i}+ e_{ij} \] + as a simple example

\[ y_{ij} = \mu + u_{j}+\beta x_{i}+ e_{ij} \]

Basic Random Effect Models

\[ y_{ij} = \mu + u_{j} + e_{ij} \]

\[ y_{ij} = \mu + u_{j} + e_{ij} \]

\[u_j \sim N(0, \sigma^2)\]

Choosing…

Forms of the random effect model

Random Intercept Model

Random Intercepts

Random Intercepts

Variance components

\(\sigma^2\) = \(\sigma^2 _{e}+\sigma^2 _{u}\) + These are called the variance components of the model, + and separate the variance into differences between individuals and differences between groups

\[\rho(y_{ij},y_{i'j}) = \frac{ \sigma^2 _{u} }{ \sigma^2 _{u} + \sigma^2 _{e} }\]

Empirical Example

In this example, I introduce how to fit the multi-level model using the lme4 package. This example considers the linear case of the model, where the outcome is assumed to be continuous, and the model error term is assumed to be Gaussian. Subsequent examples will highlight the Generalized Linear Mixed Model (GLMM).

This example shows how to: * Examine variation between groups using fixed effects * Fit the basic random intercept model :

\(y_{ij} = \mu + u_{j} + e_{ij}\) with \(u_j \sim N(0, \sigma^2)\)

Where the intercepts (\(u_j\)) for each group vary randomly around the overall mean (\(\mu\))

*I also illustrate how to include group-level covariates and how to fit the random slopes model and the model for cross level interactions

The example merges data from the 2014 CDC Behavioral Risk Factor Surveillance System (BRFSS) SMART MSA data. Link and the 2010 American Community Survey 5-year estimates at the MSA level. More details on these data are found below.

#load brfss
library(car)
library(stargazer)
library(survey)
library(sjPlot)
library(ggplot2)
library(pander)
library(knitr)
load("~/Google Drive/classes/dem7283/class18/data/brfss16_mmsa.Rdata")
set.seed(12345)
#samps<-sample(1:nrow(brfss16m), size = 40000, 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)

Recode variables

#sex
brfss16m$male<-ifelse(brfss16m$sex==1, 1, 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")

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')

#employloyment
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")

brfss16m$obese<-ifelse(is.na(brfss16m$bmi)==T, NA, 
                       ifelse(brfss16m$bmi>30,1,0))

I want to see how many people we have in each MSA in the data:

#Now we will begin fitting the multilevel regression model with the msa
#that the person lives in being the higher level
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 12260
#people within each msa

#How many total MSAs are in the data?
length(table(brfss16m$mmsa))
## [1] 143
#MSAs

Higher level predictors

We will often be interested in factors at both the individual AND contextual levels. To illustrate this, I will use data from the American Community Survey measured at the MSA level. Specifically, I use the DP3 table, which provides economic characteristics of places, from the 2010 5 year ACS Link.

library(acs)
#Get 2010 ACS median household incomes for tracts in Texas
msaacs<-geo.make(msa="*")

#ACS tables B17 is poverty, B19 is Gini, B25001 is housing vacancy, B25035 is median year built 
acsecon<-acs.fetch(key=mykey, endyear=2010, span=5, geography=msaacs, variable = c("B19083_001","B17001_001","B17001_002" , "B25002_001","B25002_003", "B25035_001"))

colnames(acsecon@estimate)
## [1] "B19083_001" "B17001_001" "B17001_002" "B25002_001" "B25002_003"
## [6] "B25035_001"
msaecon<-data.frame(gini=acsecon@estimate[, "B19083_001"], 
ppoverty=acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"],
giniz=scale(acsecon@estimate[, "B19083_001"]),
pvacant=acsecon@estimate[,"B25002_003"]/acsecon@estimate[, "B25002_001"],
ppovertyz=scale(acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"]), 
pvacantz=scale(acsecon@estimate[,"B25002_003"]/acsecon@estimate[, "B25002_001"]), 
medhouse=acsecon@estimate[, "B25035_001" ],
medhousez=scale(acsecon@estimate[, "B25035_001" ]))

msaecon$ids<-paste(acsecon@geography$metropolitanstatisticalareamicropolitanstatisticalarea)
head(msaecon)
##                                                  gini  ppoverty    giniz
## Adjuntas, PR Micro Area                         0.525 0.5925656 2.718785
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.529 0.5577311 2.847562
## Coamo, PR Micro Area                            0.532 0.5633984 2.944145
## Fajardo, PR Metro Area                          0.483 0.4309070 1.366623
## Guayama, PR Metro Area                          0.518 0.4980518 2.493425
## Jayuya, PR Micro Area                           0.502 0.5446137 1.978315
##                                                   pvacant ppovertyz
## Adjuntas, PR Micro Area                         0.1996029  6.271934
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.2005956  5.762039
## Coamo, PR Micro Area                            0.1500353  5.844994
## Fajardo, PR Metro Area                          0.3143172  3.905632
## Guayama, PR Metro Area                          0.1951184  4.888474
## Jayuya, PR Micro Area                           0.1562109  5.570031
##                                                  pvacantz medhouse
## Adjuntas, PR Micro Area                         0.9053869     1978
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.9199516     1979
## Coamo, PR Micro Area                            0.1781186     1979
## Fajardo, PR Metro Area                          2.5885036     1977
## Guayama, PR Metro Area                          0.8395882     1978
## Jayuya, PR Micro Area                           0.2687281     1978
##                                                 medhousez   ids
## Adjuntas, PR Micro Area                         0.4923898 10260
## Aguadilla-Isabela-San Sebastián, PR Metro Area 0.5982027 10380
## Coamo, PR Micro Area                            0.5982027 17620
## Fajardo, PR Metro Area                          0.3865770 21940
## Guayama, PR Metro Area                          0.4923898 25020
## Jayuya, PR Micro Area                           0.4923898 27580
summary(msaecon)
##       gini           ppoverty           giniz             pvacant       
##  Min.   :0.3510   Min.   :0.02425   Min.   :-2.88303   Min.   :0.04246  
##  1st Qu.:0.4190   1st Qu.:0.12308   1st Qu.:-0.69381   1st Qu.:0.09454  
##  Median :0.4390   Median :0.15172   Median :-0.04993   Median :0.12076  
##  Mean   :0.4406   Mean   :0.16409   Mean   : 0.00000   Mean   :0.13790  
##  3rd Qu.:0.4590   3rd Qu.:0.18746   3rd Qu.: 0.59396   3rd Qu.:0.15964  
##  Max.   :0.5730   Max.   :0.59257   Max.   : 4.26411   Max.   :0.63877  
##    ppovertyz          pvacantz          medhouse      medhousez      
##  Min.   :-2.0468   Min.   :-1.4002   Min.   :1940   Min.   :-3.5285  
##  1st Qu.:-0.6002   1st Qu.:-0.6361   1st Qu.:1968   1st Qu.:-0.5657  
##  Median :-0.1810   Median :-0.2514   Median :1975   Median : 0.1750  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   :1973   Mean   : 0.0000  
##  3rd Qu.: 0.3421   3rd Qu.: 0.3190   3rd Qu.:1980   3rd Qu.: 0.7040  
##  Max.   : 6.2719   Max.   : 7.3489   Max.   :1998   Max.   : 2.6086  
##      ids           
##  Length:955