HLM (the program) to R Conversion

trinker

ggplot2orBust
#1
I am taking a multilevel modeling class and reading Raudenbush and Bryk and using the accompanying HLM program for the class.



My assignments will need to be completed with HLM but I want to learn how the models work in R's lme4 syntax as well.

This weeks lab had us run 3 models: fully unconditional, partially conditional and fully conditional. Here is the HLM output for each:

fully unconditional
http://htmlpreview.github.io/?https://github.com/trinker/HLM2R/blob/master/fully unconditional.html

partially conditional
http://htmlpreview.github.io/?https.../HLM2R/blob/master/partially conditional.html

fully conditional
http://htmlpreview.github.io/?https://github.com/trinker/HLM2R/blob/master/fully conditional.html

I have tried to replicate this with the following R script (this will load the data and group center the SES variable):

Code Book
id: school ID number
minority: indicator of student minority race/ethnicity; 1 = minority, 0 = not minority
female: indicator of student gender; 1 = female, 0 = male
ses: standardized composite representing socioeconomic status; constructed from
variables measuring parent education, income, and occupation
mathach: mathematics achievement test score
size: total school enrollment
sector: type of school; 1 = Catholic, 0 = public
pracad: proportion of students on the academic track
disclim: a scale measuring disciplinary climate
himinty: level of minority enrollment; 1 ≥ 40%, 0 = < 40%
meanses: mean SES value for the students who are in the level-1 file and attend this school

Code:
## Load data and lme4
load(url("http://dl.dropboxusercontent.com/u/61803503/data/HSB.RData"))
head(dat)
library(lme4)

## Group Center Variables
dat$cent_ses <- with(dat, ses - ave(ses, id))

## Fully Unconditional
fu_mod <- lmer(mathach ~ 1 + (1 | id), data=dat)
summary(fu_mod)


## Partially Conditional
pc_mod <- lmer(mathach ~ cent_ses + (1 | id), data=dat)
summary(pc_mod)

## Fully Conditional
fc_mod <- lmer(mathach ~ cent_ses + (1 | id) + (1 | meanses) +(id | sector), data=dat)
summary(fc_mod)
Questions:
  1. I think I have replicated the fully unconditional and partially unconditional models correctly. How can I specify the fully conditional model in R that matches the HLM otput?
  2. How can I get the reliabilities from lme4?
  3. They use a Chi Squared test and log liklihood to look at final variance effect. How can I get lme4 to get this output?

For the Fully Conditional Model
Code:
Random Effect	Stan.Deviation	Variance	 d.f.           χ2      p-value
INTRCPT1, u0	 1.54271	 2.37996	 157	 605.29503	<0.001
SES slope, u1	 0.38590	 0.14892	 157	 162.30867	0.369
level-1, r	 6.05831	 36.70313
I assume the anova function but am unsure of which models to compare to which. I think maybe it's compared to the previous but am unsure and anova(fu_mod, pc_mod) doesn't give me the table from http://htmlpreview.github.io/?https.../HLM2R/blob/master/partially conditional.html

Thank you in advance.

My TA sent me the lab (not graded) so I could compare: http://htmlpreview.github.io/?https...r/HLM,lab1,output, including all 3 models.htm
 

trinker

ggplot2orBust
#2
Lazar gave an idea for the fully conditional mode (seen at bottom of link page), below but this doesn't match. I need fixed effects for the slopes and fixed effects for the intercepts but don't understand the syntax to get that.

Code:
load(url("http://dl.dropboxusercontent.com/u/61803503/data/HSB.RData"))
library(lme4)

## Group Center Variables
dat$cent_ses <- with(dat, ses - ave(ses, id))

lazar_mod <- lmer(mathach ~ ses + meanses + sector + (1+ses | id), data=dat)
summary(lazar_mod)

Output

Code:
> summary(lazar_mod)
Linear mixed model fit by REML ['lmerMod']
Formula: mathach ~ ses + meanses + sector + (1 + ses | id) 
   Data: dat 

REML criterion at convergence: 46547.47 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 id       (Intercept)  2.432   1.5595       
          ses          0.473   0.6877   0.27
 Residual             36.784   6.0649       
Number of obs: 7185, groups: id, 160

Fixed effects:
               Estimate Std. Error t value
(Intercept)     12.0245     0.2011   59.80
ses              2.1961     0.1224   17.94
meanses          3.1552     0.3881    8.13
sectorCatholic   1.4083     0.3091    4.56

Correlation of Fixed Effects:
            (Intr) ses    meanss
ses          0.074              
meanses      0.243 -0.260       
sectorCthlc -0.693  0.002 -0.334

I need to get 2 fixed effects for the level 2 predictors (sector and mean_ses) after regressing on group centered ses on math achievement . I know this isn't the preffered language but right now that's what I know (level 1 and level 2) and I haven't been able to move my thinking to a one equation model which lme4 seems to be based on.
 

Lazar

Phineas Packard
#3
Code:
Linear mixed model fit by REML ['lmerMod']
Formula: mathach ~ cent_ses + meanses + sector + cent_ses * meanses +      cent_ses * sector + (1 + cent_ses | id) 
   Data: dat 

REML criterion at convergence: 46503.66 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 id       (Intercept)  2.3795  1.5426       
          cent_ses     0.1013  0.3183   0.39
 Residual             36.7212  6.0598       
Number of obs: 7185, groups: id, 160

Fixed effects:
                        Estimate Std. Error t value
(Intercept)              12.0960     0.1987   60.87
cent_ses                  2.9388     0.1551   18.95
meanses                   5.3329     0.3692   14.45
sectorCatholic            1.2265     0.3063    4.00
cent_ses:meanses          1.0389     0.2989    3.48
cent_ses:sectorCatholic  -1.6426     0.2398   -6.85

Correlation of Fixed Effects:
            (Intr) cnt_ss meanss sctrCt cnt_s:
cent_ses     0.075                            
meanses      0.245  0.019                     
sectorCthlc -0.697 -0.053 -0.356              
cnt_ss:mnss  0.018  0.282  0.074 -0.026       
cnt_ss:sctC -0.052 -0.694 -0.027  0.077 -0.351
Gives same result. Note you should be using cent_ses not ses to get same results.
 

trinker

ggplot2orBust
#4
OK Lazr said throw in the interactions. Yepper that answers question 1:

Code:
load(url("http://dl.dropboxusercontent.com/u/61803503/data/HSB.RData"))
library(lme4)

## Group Center Variables
dat$cent_ses <- with(dat, ses - ave(ses, id))
lazar_mod <- lmer(mathach ~ cent_ses*meanses + cent_ses*sector + (1+ses | id), data=dat)
summary(lazar_mod)
 

spunky

Doesn't actually exist
#6
How can I get the reliabilities from lme4?
this is, unfortunately, not particularly straight-forward to do on lme4. please, keep in mind that many of the indices generated on HLM are geared towards social scientists and lme4 wasn't designed with us in mind :(

looking around Raudenbush's notes (or whatever the F his name is) i was able to find out that he defines reliability as:

\(\lambda_{j}=\frac{\tau_{00}}{\tau_{00}+\widehat{\sigma}^{2}/n_{j}}\)

which means that for every individual cluster you can calculate a reliability index \(\lambda_{i}\) the number that gets printed under Reliability estimate is the average of all the lambdas.

i thought maybe some kind soul out there may have developed a function to do this already, but the only thing i was able to find is the GmeanRel() function in the 'multilevel' package, which adds quite a few statistics that can be calculated from the output of a linear mixed effects model (such as intra-class correlations and whatnot). very fancy package, i must say... but there is one major drawback: it assumes models are fitted using nlme instead of lme4.

for example, in order to get the same reliability estimate to match your HLM output on the unconditional model you would do:

Code:
library(nlme)
library(multilevel)

tmod <- lme(mathach ~ 1, random=~1|id, data=dat)

mean(GmeanRel(tmod)$MeanRel)

> mean(GmeanRel(tmod)$MeanRel)
[1] 0.9013773
which is, i believe, what you are looking for.

now, although it seems to give you the answer you'd like for the reliability index of the partially-conditional model

Code:
library(multilevel)
library(nlme)

tmod1 <- lme(mathach ~ 1 + cent_ses, random=~1+cent_ses|id, data=dat) 
mean(GmeanRel(tmod1)$MeanRel)

[1] 0.9075635
i don't think the GmeanRel() function works for anything outside from the reliability index associated with the intercepts. you'll probably need to code your own function that fits each within-cluster regression, save the coefficients, calculate their variance and obtain the other lambda indices that HLM gives you.
 

spunky

Doesn't actually exist
#8
oh btw the answer to question three is anova(M1,M2,M3) like you suggested
but how would he test for the significance of the random effects on the intercept-only model?

sure, anova() works once you have random effects present but to test whether there is a random effect for the intercept only it will spit profanities at you. he needs to start testing from a model that has no random effects for anything to one that has one for the intercept.
 

spunky

Doesn't actually exist
#9
They use a Chi Squared test and log liklihood to look at final variance effect. How can I get lme4 to get this output?
just to be sure that you didn't miss my post in the chatbox, lme4 will not give you the chi-square test that the variance component \(H_{0}:\sigma^{2}=0\). lemme quote Faraway to explain why:

"In most cases, a test of random effects will involve a hypothesis of the form \(H_{0}:\sigma^{2}=0\). The standard derivation of the asymptotic χ2 distribution for the likelihood ratio statistic depends on the null hypothesis lying in the interior of the parameter space. This assumption is broken when we test if a variance is zero. The null distribution in these circumstances is unknown in general and we must resort to numerical methods if we wish for precise testing. If you do use the χ2 distribution with the usual degrees of freedom, then the test will tend to be conservative—the p-values will tend to be larger than they should be. This means that if you observe a significant effect using the χ2 approximation, you can be fairly confident that it is actually significant.

now...testing (nested) models where you have random effects is relatively easy with the anova() function, right? you just do anova(smaller_model, larger_model) where you only added one random effect and, based on your p-value, you can asses its significance. the first step is the tricky one because you're basically testing a model with no random effects (so a straight regression) VS one that does have random effects.

if you try to fit a model with no random effects in lme4 you get:

Code:
fu_mod <- lmer(mathach ~ 1, data=dat)
Error in mkReTrms(findbars(formula[[3]]), fr) : 
No random effects terms specified in formula
so we need to do some statistical trickery here if you wish to do everything on lme4. i'll take you by the hand and will list this in steps.

step #1

we need to fit a model (let's call it a null_model) that has no random effects and just fixed effects. in this case, that means fitting an OLS regression to your data which you will do using the lm() function. lmer() will spit profanities at you if you don't specify random effects so we'll use lm() here. this is going to be your baseline model against which you'll compare the model with a random effect for the intercept:

Code:
null_model  <- lm(mathach ~ 1, data=dat)
step #2

we're just going to fit the alternative model this time, the one with random effects and i'm gonna call it mixed_model to differentiate it

Code:
mixed_model <- lmer(mathach ~ 1 + (1 | id), data=dat)
as you can see, the difference between null_model and mixed_model is that one has a random effect for the intercept and the other one does not.

step #3

this is the tricky part because it assumes you know how a likelihood ratio test works. anova() will not work here because of the difference between fitting a model with lm() and lmer():

Code:
mixed_model <- lmer(mathach ~ 1 + (1 | id), data=dat)
null_model  <- lm(mathach ~ 1, data=dat)

anova(mixed_model, null_model)
Error in UseMethod("refitML") : 
  no applicable method for 'refitML' applied to an object of class "lm"
but that doesn't mean we cannot get the components that make up this test, i.e. the loglikelihoods. so we're gonna extract the log-likelihood for the null_model and the mixed_model as follows:

Code:
log_likelihood_null_model <- logLik(null_model)

log_likelihood_mixed_model <- logLik(mixed_model)
and we will combine them as follows (for reasons that you don't need to concern yourself with at the present moment):

Code:
diff_lik <- as.numeric(2*(log_likelihood_mixed_model-log_likelihood_null_model))

step #4

awesome! now you have the difference of log-likelihoods! and that distributes as a chi-square! so all we need to do is go and get a p-value using the chi-square as the reference distribution and that should be it!

Code:
pchisq(diff_lik,1,lower=FALSE)
[1] 9.179927e-216
which is a REALLY small number and matches the p<.001 on the HLM output, so you can conclude that the random effect for the intercept is, indeed, significant.

once you've moved on from this step, you can just use anova() to test the rest of your effects. now, we are technically still "wrong" in using the chi-square distribution as our reference distribution here, but, as Faraway said, since we got significance we can run away with it.

if you want a more precise test you can ask Jake for the bootstrapped likelihood ratio test thingy i once sent to him many-a-nights ago. he should have it somewhere.
 

Lazar

Phineas Packard
#10
Code:
m1 <- lmer(mathach ~ 1, data=dat)
m2 <- lmer(mathach ~ 1 + (1|id), data=dat)
AIC(m1); AIC(m2);
:) Not sure would have to think about it

EDIT: Just saw your post :)

EDIT 2: m1 lm not lmer
 

spunky

Doesn't actually exist
#11
Code:
m1 <- lmer(mathach ~ 1, data=dat)
m2 <- lmer(mathach ~ 1 + (1|id), data=dat)
AIC(m1); AIC(m2);
:) Not sure would have to think about it

if you fit a model without random effects in lme4 it spits profanities back at you. take a look:

Code:
> m1 <- lmer(mathach ~ 1, data=dat)
Error in mkReTrms(findbars(formula[[3]]), fr) : 
  No random effects terms specified in formula
 

spunky

Doesn't actually exist
#13
now i can go to bed... it's past 2:30am here! :p

question #2 is still somewhat open for those reliability indices for the regression coefficients. more research is needed here...
 

trinker

ggplot2orBust
#14
Thank you both for spending a great deal of time to help me match the HLM output and explaining why some aspects of the output aren't exactly best practice.

PS it's funny because I tried the anova thing using the lm model without random effects (I knew this much) but as you point out this didn't work. And I knew from Gelman and Hill (I think I read or maybe a lme4 tutorial) that you can't fit a model in lme4 w/o random effects. I didn't think about extracting the logLike, though I knew you could do this from when I studied logistic regression it had been so long I forgot you could do that.

Spunky walked us all through the night I was having and explained a typical student's stupidity.

So if I get this right testing for
with anova function is a no go because you break the assumption that "the standard derivation of the asymptotic χ2 distribution for the likelihood ratio statistic depends on the null hypothesis lying in the interior of the parameter space". The way to test this if I want to be a good boy is to use bootstrapping magic.

Also on the reliabilities where
and then you sum and divide by J, [tex]\frac{\sum\lambda_j}{J}[/tex], I was a bit surprised this wasn't more readily available in R. In meta-analysis (I may be wrong as I'm trying to connect to concepts but haven't full overlapped the two) this is a very important stat that is called Q that gets spit out right away in the metafor package as it's extremely important.

Again this thinking is most likely inaccurate, incomplete or just plain wrong but I'm processing through it all. I'm pleased I have some very knowledgeable peeps to help me work through it.
 

spunky

Doesn't actually exist
#15
So if I get this right testing for
with anova function is a no go because you break the assumption that "the standard derivation of the asymptotic χ2 distribution for the likelihood ratio statistic depends on the null hypothesis lying in the interior of the parameter space". The way to test this if I want to be a good boy is to use bootstrapping magic.
well... I wouldn't go as far as saying it's a "no go". because it has less power than a regular test you can maybe try it first. if significant, then carry on with the analysis. if non-significant (and the random effect is "large"), then proceed with the bootstrap.

In meta-analysis (I may be wrong as I'm trying to connect to concepts but haven't full overlapped the two) this is a very important stat that is called Q that gets spit out right away in the metafor package as it's extremely important.

Again this thinking is most likely inaccurate, incomplete or just plain wrong but I'm processing through it all. I'm pleased I have some very knowledgeable peeps to help me work through it.
have you tried to fit your models using the meta package? I don't do metanalysis, but you may be right. maybe you just need to re-frame these models in some way the metafor package can estimate them so you get that Q statistic.
 

trinker

ggplot2orBust
#16
have you tried to fit your models using the meta package? I don't do metanalysis, but you may be right. maybe you just need to re-frame these models in some way the metafor package can estimate them so you get that Q statistic.
Hmm interesting idea.
 

trinker

ggplot2orBust
#17
@Spunky why do you multiple the differences between models by 2 in this post:

Code:
diff_lik <- as.numeric(2*(log_likelihood_mixed_model-log_likelihood_null_model))
 

trinker

ggplot2orBust
#20
Chatbox info on p-values that Doug Bates detests:

trinker said:
So if one knew it was wrong to calculate pvalues but still did and the had N = 8000, J - 348...
Jake said:
but a rough rule of thumb you could use for a 2-level hierarchical model is that if the predictor varies within schools, N - J will be close, while if it varies within schools, J will be close (but note both of these ignore number of parameters)
spunky said:
JAke! we summon thee and thy power to correct our misguided p-values!
Jake said:
you can use lmerTest for welch-satterthwaite, or pbkrtest for kenward-roger. i prefer the latter package but both are good
there are also approaches not based on degrees of freedom;
see the R glmm FAQ wiki;
spunky said:
non-degrees-of-freedom approach to significance testing? what is this sorcery you speak of?
Jake said:
note both of these are only valid for linear mixed models (not generalized). most of the other approaches work on both LMMs and GLMMs
Dason said:
@spunky - permutation/randomization tests don't rely on degrees of freedom
Jake said:
all kinds of stuff spunky. parametric bootstrap, LRT, MCMC (your favorite?!), Wald test, ...