# Multilevel regression model using “multilevel” and “lme4” R packages?

After all data were collected, statistical analyses were performed to test for correlations between TI and D_d, RI, RV, and MRP. Due to the high amount of uncertainty introduced when comparing TI estimates between Groups, typical regression analyses cannot be conducted on these data. A typical regression model requires that the data are independent of each other. However, the collected data within a group are dependent on the error associated with that individual group. Therefore, a typical regression model would produce inaccurate results. Consequently, I have to use a set of multilevel regression models, which is appropriate for nested data.

In this case, estimated TI values are nested within Group, and using a multilevel model allowed us to compare TI both within and between scenes. This comparison was viable even when the data exhibit different slopes and y-intercepts caused by variations in uncertainty between group. I have to use the multilevel and lme4 packages with the R statistical language to run multilevel regression models on all sets of data. After the models were run, the statistical significance of each model result was determined using a likelihood ratio test to attain p-values.

Since I am new to R and statistics, if anyone could help me in understanding the design it will be of great help to me. I am attaching a screenshot of my dataframe. I really need help. ID = 138 Group = 18

I have to account for nesting by both random intercept and random slopes, for a single crater ID, multiple values, i.e., TI,RI,RV,D_d and MRP is measured. example: for ID 103, TI, RI, RV, D_d and MRP is measured, similarly for each crater these parameters were measured. From what I can see of your data and the descriptions, you do not have multiple measures within ID. You have measured several variables, D_d, RI, RV, and MRP once for each ID.

Thus ID seems is the unit of measurements (that is, it unique to each row in your data).

However you do seem to have multiple measures within Group, and therefore a model with random intercepts for Group would seem to be appropriate. I would therefore suggest the following model as a starting point:

lmer(TI ~ D_d + RI + RV + MRP + (1 | Group), data = ... )


This will estimate fixed effects for D_d, RI, RV, and MRP, along with a variance for the random Group variable, which will account for the non-independence of measurements within each group.

• Here ID represents the crater number, the craters are organized by the Group that cover them. example: 5 craters cover group A. I tried the model described by you, but I am not getting the correct p-values. I don't know where I am going wrong. – Farzana Feb 4 at 13:25
• @Farzana OK that is what I thought, so ID is nested in Group, but each ID/crater has it's measurements taken only once. That is, each row represents one crater. – Robert Long Feb 4 at 13:26
• Each ID consists of measurement of 5 cases (i.e. TI, RV,RI,D_d and MRP), error were determined based on repeat measurements of a particular crater. – Farzana Feb 4 at 13:32
• @Farzana so you measured TI, RV,RI,D_d and MRP 5 times for each crater ? Or you measured them all once for each crater – Robert Long Feb 4 at 13:38
• No, I measured it (TI,RI,RV,D_d,MRP) once for each crater. – Farzana Feb 5 at 3:08

A couple of points:

• Mixed models are indeed used to account for correlations in your outcome variable, I guess TI within the levels of grouping/cluster variables, i.e., ID and Group in your case. Assuming that normal error terms would be adequate for TI, you could use a linear mixed model. For example, using function lmer() from package lme4, e.g.,

fm1 <- lmer(TI ~ RI + (1 | Group / ID), data = tisia)


If you are going also to load the lmerTest package, you will obtain a p-value for the association between TI and RI.

• Model fm1 above postulates that the correlation between any pair of measurements of TI within the same combination of levels of ID and Group is the same. If instead, you want to assume that the correlations within the same combination of ID and Group decay as the difference in RI values increases, then you could include the random slope for RI, i.e.,

fm2 <- lmer(TI ~ RI + (RI | Group / ID), data = tisia)

• You could compare the two models to see if this improves the fit using a likelihood ratio test implemented by the anova() function, i.e.,

anova(fm1, fm2)

• Thank you for your response, I followed your instruction and tried applying those models in r, but an error message pops up which states : Error: number of levels of each grouping factor must be < number of observations, how to resolve this error? – Farzana Feb 4 at 12:38
• Couple of things: 1) Did you code your two grouping variables in such a way that they are distinct? Imagine that in Group A, there are 5 IDs, 100-104, and then in Group B, there also 5 IDs, 100-104. This shouldn't be a problem if you have the nesting correct, but it's always a good idea to give IDs unique and distinct values across groups. 2) Try running the random effects as (1|Group) + (1|ID). What happens? Are the reported Ns of each of these categories correct? – Erik Ruzek Feb 4 at 21:05