10  Multilevel plm analyses

Note

Read Chapter 9 before you start with this chapter.

The hplm function call

hplm(data, dvar, pvar, mvar, model = c(“W”, “H-M”, “B&L-B”, “JW”), contrast = c(“first”, “preceding”), contrast_level = NA, contrast_slope = NA, method = c(“ML”, “REML”), control = list(opt = “optim”), random.slopes = FALSE, lr.test = FALSE, ICC = TRUE, trend = TRUE, level = TRUE, slope = TRUE, fixed = NULL, random = NULL, update.fixed = NULL, data.l2 = NULL, …)

Multilevel analyses can take the piecewise-regression approach even further. It allows for

The basic function for applying a multilevel piecewise regression analysis is hplm. The hplm function is similar to the plm function, so I recommend that you get familar with plm before applying an hplm.

Here is a simple example:

hplm(exampleAB_50)
Hierarchical Piecewise Linear Regression

Estimation method ML 
Contrast model: W / level: first, slope: first
50 Cases

ICC = 0.287; L = 339.0; p = 0.000

Fixed effects (values ~ 1 + mt + phaseB + interB)

                   B    SE   df      t p
Intercept     48.398 1.484 1328 32.611 0
Trend mt       0.579 0.116 1328  5.006 0
Level phase B 14.038 0.655 1328 21.436 0
Slope phase B  0.902 0.119 1328  7.588 0

Random effects (~1 | case)

          EstimateSD
Intercept      9.970
Residual       5.285

Here is an example inlcuding random slopes:

hplm(exampleAB_50, random.slopes = TRUE)
Hierarchical Piecewise Linear Regression

Estimation method ML 
Contrast model: W / level: first, slope: first
50 Cases

ICC = 0.287; L = 339.0; p = 0.000

Fixed effects (values ~ 1 + mt + phaseB + interB)

                   B    SE   df      t p
Intercept     48.211 1.398 1328 34.497 0
Trend mt       0.621 0.113 1328  5.516 0
Level phase B 13.872 0.894 1328 15.513 0
Slope phase B  0.864 0.116 1328  7.433 0

Random effects (~1 + mt + phaseB + interB | case)

              EstimateSD
Intercept          9.352
Trend mt           0.096
Level phase B      4.537
Slope phase B      0.126
Residual           4.974

10.0.1 Adding additional L2-variables

The add_l2 function call

add_l2(scdf, data_l2, cvar = “case”)

In some analyses researchers want to investigate whether attributes of the individuals contribute to the effectiveness of an intervention. For example might an intervention on mathematical abilities be less effective for student with a migration background due to too much language related material within the training. Such analyses can also be conducted with scan. Therefore, we need to define a new data frame including the relevant information of the subjects of the single-case studies we want to analyze. This data frame consists of a variable labeled case which has to correspond to the case names of the scfd and further variables with attributes of the subjects. To build a data frame we can use the R function data.frame.

L2 <- data.frame(
  case = c("Antonia","Theresa", "Charlotte", "Luis", "Bennett", "Marie"), 
  age = c(16, 13, 13, 10, 5, 14), 
  sex = c("f","f","f","m","m","f")
)
L2
       case age sex
1   Antonia  16   f
2   Theresa  13   f
3 Charlotte  13   f
4      Luis  10   m
5   Bennett   5   m
6     Marie  14   f

Multilevel analyses require a high number of Level 2 units. The exact number depends on the complexity of the analyses, the size of the effects, the number of level 1 units, and the variability of the residuals. But surely we need at least about 30 level 2 units. In a single-case design that is, we need at least 30 single-cases (subjects) within the study. After setting the level 2 data frame we can merge it to the scdf with the add_l2() function (alternatively, we can use the data.l2 argument of the hplm function). Then we have to specify the regression function using the update.fixed argument. The level 2 variables can be added just like any other additional variable. For example, we have added a level 2 data-set with the two variables sex and age. update could be construed of the level 1 piecewise regression model .~. plus the additional level 2 variables of interest + sex + age. The complete argument is update.fixed = .~. + sex + age. This analyses will estimate a main effect of sex and age on the overall performance. In case we want to analyze an interaction between the intervention effects and for example the sex of the subject we have to add an additional interaction term (a cross-level interaction). An interaction is defined with a colon. So sex:phase indicates an interaction of sex and the level effect in the single case study. The complete formula now is update.fixed = .~. + sex + age + sex:phase.

scan includes an example single-case study with 50 subjects example50 and an additional level 2 data-set example50.l2. Here are the first 10 cases of example50.l2.

case sex age
Roman m 12
Brennen m 10
Ismael m 13
Donald m 11
Ricardo m 13
Izayah m 11
Ignacio m 12
Xavier m 12
Arian m 10
Paul m 10

Analyzing the data with hplm could look like this:

exampleAB_50 %>%
  add_l2(exampleAB_50.l2) %>%
  hplm(update.fixed = .~. + sex + age)
Hierarchical Piecewise Linear Regression

Estimation method ML 
Contrast model: W / level: first, slope: first
50 Cases

ICC = 0.287; L = 339.0; p = 0.000

Fixed effects (values ~ mt + phaseB + interB + sex + age)

                   B     SE   df      t     p
Intercept     44.878 11.926 1328  3.763 0.000
Trend mt       0.581  0.116 1328  5.026 0.000
Level phase B 14.023  0.655 1328 21.405 0.000
Slope phase B  0.900  0.119 1328  7.569 0.000
sexm          -6.440  2.727   47 -2.362 0.022
age            0.603  1.073   47  0.562 0.577

Random effects (~1 | case)

          EstimateSD
Intercept      9.446
Residual       5.284
# Alternatively:
# hplm(exampleAB_50, data.l2 = exampleAB_50.l2, update.fixed = .~. + sex + age)

sex is a factor with the levels f and m. So sexm is the effect of being male on the overall performance. age does not seem to have any effect. So we drop age out of the equation and add an interaction of sex and phase to see whether the sex effect is due to a weaker impact of the intervention on males.

exampleAB_50 %>%
  add_l2(exampleAB_50.l2) %>%
  hplm(update.fixed = .~. + sex + sex:phaseB)
Hierarchical Piecewise Linear Regression

Estimation method ML 
Contrast model: W / level: first, slope: first
50 Cases

ICC = 0.287; L = 339.0; p = 0.000

Fixed effects (values ~ mt + phaseB + interB + sex + phaseB:sex)

                        B    SE   df       t    p
Intercept          48.573 1.968 1327  24.676 0.00
Trend mt            0.609 0.109 1327   5.573 0.00
Level phase B      17.726 0.684 1327  25.922 0.00
Slope phase B       0.884 0.112 1327   7.868 0.00
sexm               -0.593 2.741   48  -0.216 0.83
Level phase B:sexm -7.732 0.609 1327 -12.699 0.00

Random effects (~1 | case)

          EstimateSD
Intercept      9.494
Residual       4.989

Now the interaction phase:sexm is significant and the main effect is no longer relevant. It looks like the intervention effect is \(7.7\) points lower for male subjects. While the level-effect is \(17.7\) points for female subjects it is \(17.7\) - \(7.7\) = \(10\) for males.