# Variable selection in Hierarchical Linear Modelling HLM through nlme lme()

Background of my question:- In Linear Regression through R we can mention the direction="both"/"forward"/"backward" in step(lm()) function to tell R for choosing the best set of variables based on AIC. The output which we get is the final selection of a reduced number of attributes that best explains my dependent variable. My Question:- If I want to use HLM method through nlme package for a dataset with many attributes then is there any optimized method to choose best selection of variables? E.g code:-

model_mydf <- lme(lgaspcar ~ 1 + lincomep+lrpmg+lcarpcap,
data = Gasoline,
method = "ML",control = ctrl,
random = ~ 1 +lincomep +lrpmg + lcarpcap | country)


This data comes from data("Gasoline",package = "plm"). In this case I have only 3 explanatory variables and it is less time consuming to iterate and see which is my best set of variables. Suppose I have larger attributes (or explanatory variables) then we need to keep on iterating lme() function to check the coefficients, p-value, Rsquare and MAPE. Because there is no option to mention direction="" argument similar to step(lm()). How can I reduce my iteration time if I want to adopt HLM method through lme()? Is there anything similar to step-wise/backward/forward in lme() function?

With response to @Gregg H dated 31st March 2018:

ctrl<- lmeControl(opt='optim',optimMethod="SANN")
set.seed(10)
fe_ML <- lme(lgaspcar ~ 1,
data = Gasoline,
method = "ML",control = ctrl,
random = ~1 | country)
summary(fe_ML)
D_ML <- 2*fe_ML$logLik df_ML <- fe_ML$fixDF$terms[] ctrl<- lmeControl(opt='optim',optimMethod="SANN") set.seed(10) re_REML <- lme(lgaspcar ~ 1, data = Gasoline, method = "REML",control = ctrl, random = ~1 | country) D_REML <- 2*re_REML$logLik
df_REML <- re_REML$fixDF$terms[]


D_ML & D_REML are respective 2*LL for fixed effect and random effect models using ML and REML method. df_ML & df_REML are respective degrees of freedom (in this case both are 324). How do I do a chisquare test now using chisq.test()? How do I keep iterating for lincomep+lrpmg+lcarpcap one by one to get the best combination? In above case I have tested only intercept. How do I combine both ML and REML combination of variables in one final lme() function to get the best combination result?

With response to @BenBolker dated 31st March 2018:
ctrl<- lmeControl(opt='optim',optimMethod="SANN")
set.seed(10)
f <- lgaspcar ~ 1
r <- ~ 1
model_mydf <- lme(f,
data = mydf,
method = "ML",control = ctrl,
random = r | country)

• Please clarify what exactly you want to change in the model selection process: the fixed effects or the random effects. You need different approaches for each. – Gregg H Mar 30 '18 at 18:08
• Hi @GreggH , Thank you for answering. I want to check both. I understand there are some attributes which would not vary across entity/time component in HLM and they would remain in fixed effect and not in random effect. I also understand that in lme() function if I keep increasing the variables in random= argument then processing time increases manifold and we need to avoid that. However there would be some attributes in the data which would remain in fixed and random effects both. In those cases how should we approach best variable selection process in lme() alike step(lm())? Please HELP! – ARIMITRA MAITI Mar 30 '18 at 19:01
• if you fit with ML you can use drop1() to test one-at-a-time deletion ... – Ben Bolker Mar 30 '18 at 21:30
• HI @BenBolker. I tried as per suggestion with add1() instead and I do get results similar to stepAIC() from MASS library. Please find my try on "With response to BenBolker dated 31st March 2018:" in the original question. How do I keep adding attributes now based on AIC values. The random = argument is not accepting any objects unless I put the attributes manually in there. Please help. – ARIMITRA MAITI Apr 1 '18 at 6:19

First, I do not know if there is a comparable set of packages or functions like can do automated model comparison testing (similar to step() or regsubsets()).

Second, if you are comparing changes in the random effects, you need to use REML.  If you are comparing changes in fixed effects, you need to use ML.

With those observations out of the way, you can compare two nested models using the deviance score, which is $-2$ times the log-likelihood: $D = -2·LL$.  $\Delta D = D_s - D_b$ (where $s$ indicates smaller model and $b$ indicates the bigger model).  This value should follow a chi-squared distribution where the degrees of freedom are the difference in number of parameters in the models.  If $∆D$ is statistically significant, then the models are different from each other.

Hope this helps.

As requested, here is some updated code (revised from that provided in the question) to demonstrate the chi-square test for these type of model comparisons:

library(nlme)
data("Gasoline",package = "plm")

ctrl <- lmeControl(opt='optim',optimMethod="SANN")
set.seed(10)
fe_ML <- lme(lgaspcar ~ 1,
data = Gasoline,
method = "ML",control = ctrl,
random = ~1 | country)
summary(fe_ML)
D_ML <- -2*fe_ML$logLik df_ML &lt;- fe_ML$fixDF$terms[] fe_ML2 <- lme(lgaspcar ~ 1 + I(year-1960), data = Gasoline, method = "ML",control = ctrl, random = ~1 | country) summary(fe_ML2) D_ML2 <- -2*fe_ML2$logLik
df_ML2 &lt;- fe_ML2$fixDF$terms[]

cat("\nchi-square statistic:  ");  cat(D_ML-D_ML2)
cat("\ndegres of freedom:  ");  cat(df_ML-df_ML2)
cat("\nP-value:  ");   cat(pchisq(D_ML-D_ML2,df_ML-df_ML2,lower.tail=FALSE))


• Would you be generous to render some more help on my further queries? I have pasted in original question "With response to @GreggH dated 31st March 2018" Please help! – ARIMITRA MAITI Apr 1 '18 at 6:06
• I will add code to my answer (Addendum #1) – Gregg H Apr 1 '18 at 12:56
• I do not know if there is an update() function comparable for lme(). (There may be, so it is probably worth exploring.) A brute force approach is to create a string variable with the formula, and add "+new.var" to the string, and then use the as.formula() function. – Gregg H Apr 1 '18 at 13:03
• Thank you ever so much. I will come back soon to seek your review on my progress. – ARIMITRA MAITI Apr 1 '18 at 13:54