0
$\begingroup$

I've built a lmertree model using the GLMERTREE package with random slope and intercept, a treatment variable, and a "partitioning" variable. I've attached a toy dataset and the associated model:

Dataset:

toy_data<- structure(list(abund_per_night = c(3.5572510116355, 4.01750851327412, 
3.46350126499831, 3.46332916276916, 3.63946409455059, 4.94499813045628, 
4.22468459335564, 4.0623568867287, 2.88240358824699, 4.53259949315326, 
3.5727776966247, 2.70805020110221, 3.11905548958599, 3.17108516103185, 
1.32175583998232, 3.01343185065339, 4.15612192191834, 3.95603989084492, 
2.55128171873287, 3.97582596198676, 2.36712361413162, 2.83321334405622, 
2.28238238567653, 1.88273124743378, 1.33500106673234, 2.70327692234955, 
3.57531544700328, 4.57814151360017, 1.87180217690159, 2.66258782702545, 
3.14127457884465, 2.0636931847117, 2.3434070875143, 2.23359222150709, 
3.86223234092513, 3.6329313801814, 3.21386328304466, 2.05713578416554, 
3.01099974568478, 1.62362254742606, 2.10413415427021, 1.73460105538811, 
3.54095932403731, 1.36524095192206, 1.70474809223843, 2.25438299117617, 
2.09714111877924, 1.56397553835734, 2.18925640768704), Site.Code = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L
), .Label = c("Site3", "Site4", "Site6", "Site7", "Site8"), class = "factor"), 
    pdsi_500_Jul = c(0.75385445356369, 3.04974412918091, 6.07198905944824, 
    4.90420579910278, 4.46486043930054, -3.06367492675781, 3.16913890838623, 
    2.77952671051025, 2.4449200630188, 3.67740941047668, 0.988260984420776, 
    2.85133194923401, 5.53220462799072, 5.27635478973389, 4.83563280105591, 
    -2.81783366203308, 2.76242399215698, 2.30299067497253, 2.59277606010437, 
    3.73904895782471, 0.869111835956573, 2.70135831832886, 6.06647205352783, 
    5.31293201446533, 4.80502414703369, -2.98165917396545, 2.97572326660156, 
    2.55066156387329, 2.53725910186768, 4.15768718719482, 0.905851364135742, 
    5.67419719696045, 5.34246158599854, 4.60490322113037, -3.07186555862427, 
    2.56848764419556, 2.48435306549072, 2.96440410614014, 4.11801624298096, 
    1.48729205131531, 2.74770283699036, 5.51829481124878, 5.98142242431641, 
    4.98692893981934, -2.68818712234497, 2.20900440216064, 3.89873313903809, 
    3.3214693069458, 3.74352169036865), forusgs_1000 = c(16.6922033898305, 
    16.6922033898305, 16.6922033898305, 16.6922033898305, 16.6922033898305, 
    16.6922033898305, 16.6922033898305, 16.6922033898305, 16.6922033898305, 
    16.6922033898305, 32.8421052631579, 32.8421052631579, 32.8421052631579, 
    32.8421052631579, 32.8421052631579, 32.8421052631579, 32.8421052631579, 
    32.8421052631579, 32.8421052631579, 32.8421052631579, 11.1906377204885, 
    11.1906377204885, 11.1906377204885, 11.1906377204885, 11.1906377204885, 
    11.1906377204885, 11.1906377204885, 11.1906377204885, 11.1906377204885, 
    11.1906377204885, 0.831030577576444, 0.831030577576444, 0.831030577576444, 
    0.831030577576444, 0.831030577576444, 0.831030577576444, 
    0.831030577576444, 0.831030577576444, 0.831030577576444, 
    7.95168972556135, 7.95168972556135, 7.95168972556135, 7.95168972556135, 
    7.95168972556135, 7.95168972556135, 7.95168972556135, 7.95168972556135, 
    7.95168972556135, 7.95168972556135)), .Names = c("abund_per_night", 
"Site.Code", "pdsi_500_Jul", "forusgs_1000"), row.names = c(21L, 
22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 
35L, 36L, 37L, 38L, 39L, 40L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 
58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 349L, 
350L, 351L, 352L, 353L, 354L, 355L, 356L, 357L, 358L), class = "data.frame")

Model:

lmer.tree<-lmertree(abund_per_night ~ pdsi_500_Jul | (1+pdsi_500_Jul|Site.Code) | forusgs_1000 , data = toy_data, joint=F)

First, I'm trying to determine how to calculate confidence intervals for the regression coefficients at each of the terminal nodes, which I can identify using the coef() function.

coef(glmer.tree)

I'm also considering a model that includes a fixed effect in the lmer/random part of the lmertree formula.

lmer.tree2<-lmertree(abund_per_night ~ pdsi_500_Jul | pdsi_500_Jul + (1+pdsi_500_Jul|Site.Code) | forusgs_1000 , data = toy_data, joint=F)

This allows me to generate confidence intervals on each of the random slopes using the below code, which I'm doing because I want to look at the distribution and uncertainty in random slope estimates to address my research questions. Although the code runs with the fixed effect in the "random" part of the formula, I'm not sure the package is designed to run in that way and I'm curious what the implications are for the other parameter estimates. I'm also curious why if many of the random slopes are positive and the vast majority of CIs on the random slopes overlap 0, the coefficient estimates at the terminal nodes are both strongly negative. Any insight into these questions would be greatly appreciated.

rslopes <- coef(lmer.tree2$lmer)$Site.Code[,2]
varfix <- vcov(lmer.tree2$lmer)[2,2]
re <- ranef(lmer.tree2,condVar=TRUE)
varcm <- attr(re$Site.Code,"postVar")[2,2,]
vartot <- varfix+varcm
plotCI(1:length(rslopes),rslopes,2*sqrt(vartot))
upper_CI<-rslopes + 2*sqrt(vartot)
lower_CI<-rslopes - 2*sqrt(vartot)
Improve this question
$\endgroup$
2
  • $\begingroup$ Although this has a lot of R in it, I think the underlying question really is statistical, not programming; I think it's more likely to find knowledgeable answerers here than on SO, so I am voting to leave it open. $\endgroup$
    – Peter Flom
    Feb 27, 2018 at 12:21
  • 1
    $\begingroup$ Agreed, thanks for doing so, @PeterFlom. I tried to focus on the statistical aspects in my answer. $\endgroup$
    – Achim Zeileis
    Feb 28, 2018 at 21:22

1 Answer 1

Reset to default
2
$\begingroup$

The glmertree package estimates (generalized) linear mixed effect models of the form: $$ \begin{eqnarray*} \mathrm{E}(y_{i}) & = & \mu_{i}\\ g(\mu_{i}) & = & x_i^\top \beta_{j(i)} + z_i^\top b \end{eqnarray*} $$ where $j(i)$ is the subgroup (or terminal node or leaf) $j$ that observation $i$ gets assigned to by the tree. (For more details see Fokkema et al. 2018, Behavior Research Methods, doi:10.3758/s13428-017-0971-x, Equation 5 where a slightly different notation is used).

Thus, it is possible to estimate the whole tree in a single "joint" GLMM once the group structure $j(i)$ is "known" (actually estimated). The estimate can be seen as an interaction model where all fixed effects with regressors $x_i$ interact with a factor coding the grouping structure whereas a global estimate is used for the random effects with regressors $z_i$.

This view is also implemented by default in glmertree. For some reason you have set joint = FALSE which we have added because other approaches in the literature do this. But it is really inferior because: (a) the EM-style iteration needs longer to converge, (2) yet convergence is not as close to the maximum, (3) you do not get a single model object containing all parameters. If you use the default joint = TRUE you get:

tr <- lmertree(abund_per_night ~ pdsi_500_Jul | (1 + pdsi_500_Jul | Site.Code) | forusgs_1000 , data = toy_data)
tr$lmer
## Linear mixed model fit by REML ['lmerMod']
## Formula: abund_per_night ~ .tree + .tree:pdsi_500_Jul + (1 + pdsi_500_Jul |  
##     Site.Code)
##    Data: data
## Weights: .weights
## REML criterion at convergence: 117.2104
## Random effects:
##  Groups    Name         Std.Dev. Corr 
##  Site.Code (Intercept)  0.6991        
##            pdsi_500_Jul 0.1005   -1.00
##  Residual               0.7022        
## Number of obs: 49, groups:  Site.Code, 5
## Fixed Effects:
##         (Intercept)               .tree3  .tree2:pdsi_500_Jul  
##             2.68429              1.09101             -0.06950  
## .tree3:pdsi_500_Jul  
##            -0.09164

In principle, it is technically possible to obtain confidence intervals from this object. However: Doing so ignores that the subgroup structure from the tree is not exogenously known in advance but actually estimated from the data. This is similar to performing inference after model selection via AIC or BIC. The latter is often done in practice but it is also known from the literature that it can go pathologically wrong (see the work by Leeb and Poetscher, e.g., their 2005 Econometric Theory paper "Model Selection and Inference: Facts and Fiction", doi:10.1017/S0266466605050036, is a good starting point). If you want to have a look nevertheless you can do so by:

confint(tr$lmer)
## Computing profile confidence intervals ...
##                           2.5 %     97.5 %
## .sig01               0.07944129 1.28451505
## .sig02              -1.00000000 1.00000000
## .sig03               0.00000000 0.22085933
## .sigma               0.56444388 0.85935443
## (Intercept)          1.83700696 3.52714904
## .tree3              -0.24063544 2.42582741
## .tree2:pdsi_500_Jul -0.22397085 0.08471932
## .tree3:pdsi_500_Jul -0.28142691 0.09911076

This throws a number of warnings, though, probably associated with the awkward results in .sig02 and .sig03. Possibly this indicates that the model is already quite complex for the size of the toy_data.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for the clarification! I understand your point regarding estimation of confidence intervals. 1.) Is the solution to subset the data and generate CIs on a test set? 2.) Or couldn’t some form of cross validation be used to create these estimates/are you aware of a package that can do this with a lmertree object? There has to be some way to generate unbiased uncertainty estimates on these values, otherwise aren't they are essentially useless? 3.) Is there a way to generate CIs on the global random slope estimates, even if they suffer from the same issues? $\endgroup$
    – Nick_89
    Feb 28, 2018 at 19:55
  • $\begingroup$ Also, just out of curiosity, is it possible that the tree algorithm creates as many subgroups j as there are i such that each terminal node is the same as the random slope/intercept? How would you interpret this scenario? $\endgroup$
    – Nick_89
    Feb 28, 2018 at 19:56
  • $\begingroup$ (1) Probably better but I'm not sure it gives valid CIs in all data-generating processes. (2) See stats.stackexchange.com/questions/328018/… (3) The CIs on the global estimates are probably much more robust but I'm not sure whether they are valid asymptotically. Note also that the fact that there is (to the best of my knowledge) no unifying asymptotic theory for trees and things can go wrong does not mean that they must go wrong. I find the CIs useful as "yardsticks" albeit I would not use them for confirmatory inference. $\endgroup$
    – Achim Zeileis
    Feb 28, 2018 at 20:15
  • $\begingroup$ And if you know the subgroups (i = j) then you don't need to estimate a tree?! $\endgroup$
    – Achim Zeileis
    Feb 28, 2018 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.