# GLMERTREE Confidence interval

I am currently analysing a small dataset (see sample data below) using lmertree.

My code:

a <-lmertree(Percent ~ 1 | AREA | A + B, data=sample_data)
plot(a, which = "tree")

The lmertree result (plot)

I have 3 questions:

1. I tried getting the confidence intervals of the split values (39.494 & 0.397) with the package "stablelearner", I cannot get functions like "stabiliy" to work (Learner of class: stabletree not in LearnerList.) , I can only get following function to work:

st <- stabletree(a)

which gives me the result:

B = 500
Method = Bootstrap sampling with 100.0% data

Variable selection overview:

freq *  mean *
B 0.678 1 0.876 1
A 0.592 1 0.650 1
(* = original tree)

I am not sure how to interpret these though, or how to calculate the 95% confidence interval for the split values.

1. Is it right to interpret the lmertree result: Predictor "A" is more important in explaining the variance in "Percent" as "B", because it is the first split of the tree?

2. A & B are correlated-Spearman rank: 0.52, is this a problem for the lmertree?

Sample data: (.csv)

AREA;Year;Percent;A;B
1,00;2014;16,374;57,310;0,385
1,00;2015;16,364;56,970;0,392
1,00;2016;13,333;56,970;0,385
1,00;2017;4,819;56,024;0,386
2,00;2013;11,765;0,000;0,345
2,00;2014;5,911;0,000;0,345
2,00;2015;10,837;0,000;0,394
2,00;2016;23,469;0,000;0,393
2,00;2017;11,330;0,000;0,394
3,00;2005;14,740;47,010;0,293
4,00;2013;4,158;55,799;0,186
4,00;2014;2,655;55,752;0,225
4,00;2015;3,571;56,818;0,145
4,00;2016;3,873;54,930;0,145
5,00;2013;16,848;27,174;0,524
5,00;2014;10,149;29,455;0,521
5,00;2015;15,591;27,151;0,577
5,00;2016;27,390;29,199;0,572
6,00;2016;7,614;30,964;0,291
6,00;2017;33,712;30,682;0,289
7,00;2013;19,091;0,000;0,397
7,00;2014;38,952;0,000;0,433
7,00;2015;37,011;0,000;0,432
8,00;2014;40,796;25,373;0,458
8,00;2015;43,781;25,373;0,458
9,00;2013;10,606;82,576;0,423
9,00;2014;10,680;77,670;0,368
9,00;2015;8,911;77,723;0,253
9,00;2016;6,633;76,531;0,260
10,00;2013;31,658;0,000;0,457
10,00;2014;48,454;0,000;0,454
10,00;2015;22,222;0,000;0,456
11,00;2014;39,747;39,494;0,554
11,00;2015;44,020;39,440;0,558
12,00;2005;9,111;50,000;0,203
12,00;2006;11,233;50,889;0,215
12,00;2007;8,333;50,889;0,207
12,00;2008;11,313;52,889;0,196
12,00;2009;10,805;52,444;0,199
12,00;2010;7,692;53,111;0,191
12,00;2011;8,300;52,889;0,177
12,00;2012;11,395;53,111;0,210
12,00;2013;5,788;51,778;0,212
12,00;2014;11,546;52,444;0,196
13,00;2012;36,000;24,500;0,383
13,00;2013;39,000;24,500;0,383
13,00;2014;47,000;24,500;0,557
13,00;2015;31,156;24,623;0,558
13,00;2016;55,330;23,858;0,555
13,00;2017;38,384;24,242;0,556
• It's a good idea to tag questions with the language you use (R ???) if only to attract or repel people here according to their expertise. In this case, Q1 looks off-topic here (getting code to work) while Q2 and Q3 do seem to have statistical flavour. I will not vote to close as off-topic as too close to a software-specific question, but it's open to people closer to this method to judge on that. – Nick Cox Feb 11 '18 at 11:44
• @NickCox I feel that the question is more appropriate for this forum as opposed to SO. While the questions are clearly connected to our R packages glmertree (on CRAN) and stablelearner (on R-Forge), the problems tackled are more of a conceptual/statistical nature. Re: CSV. A semicolon-separated file is produced on many European-language Excel systems because the comma instead of the period is used in these languages as the decimal separator. For this reason, R offers both read.csv() (sep="," and dec=".") and read.csv2() (sep=";" and dec=","). – Achim Zeileis Feb 11 '18 at 13:28
• @AchimZeileis Thanks for clarification. OP gets attention from the top authority; that's excellent. (I do understand about different separators, particularly if commas are decimal separators.) – Nick Cox Feb 11 '18 at 13:34

1. You do not get classical bootstrap confidence intervals for trees. Already in the first split not only the split point but also the split variable might vary (especially if the variables are correlated). Subsequent splits may then pertain to completely different earlier partitions and hence cannot be compared easily. Therefore, the result of stabletree() gives more guidance about the stability of various properties of the tree and not so much classical inference. See the corresponding manuscript for more details, a preprint version is available on my web page. https://eeecon.uibk.ac.at/~zeileis/papers/Philipp+Zeileis+Strobl-2016.pdf
2. If you inspect the results more closely you will see that B is indeed used more often as a splitting variable but that A is also used frequently. Furthermore, if you look at image(st) you see that there are quite a few bootstrap samples where B but not A - and also A but not B is used. So both factors appear relevant to some degree.
3. As the variables are considered one at a time in the tree, there are no numerical problems. However, as already discussed above, you see that it is hard to separate in a data-driven way which of the two variables is "more important".

Additionally, I would like to point out:

• There is quite a bit of heteroscedasticity in the response which is not surprising for something in percent. As your tree shows the variance is clearly lower for low percentages as opposed to percentages around 50%. To stabilize the variance a logit transformation is often used.
• All analyses from simple scatterplots over the GLMM tree to the corresponding stability analysis suggest that there is a rather linear relationship with B but more of a step-shaped influence of A. Maybe there is either more data or subject-matter knowledge that would support this.

Illustration: Read data, transformations, and exploratory display.

d <- transform(d,
AREA = factor(AREA),
Logit = qlogis(Percent/100)
)
plot(Logit ~ A, data = d)
plot(Logit ~ B, data = d)

These show a step at about A = 40 (left) and more of a linear relationship in B (right).

The corresponding GLMM tree yields two splits in B. Note, however, that this seems to approximate the linear pattern seen in the scatterplot above.

tr <- lmertree(Logit ~ 1 | AREA | A + B, data = d)
plot(tr, which = "tree")

The stability analysis shows that not only B is used as a splitting variable on bootstrap samples but also A. If A is used the split point is mostly around 40. If B is used there is more variable in the split points.

set.seed(1)
st <- stabletree(tr)
summary(st)
## Call:
## lmertree(formula = Logit ~ 1 | AREA | A + B, data = d)
##
## Sampler:
## B = 500
## Method = Bootstrap sampling with 100.0% data
##
## Variable selection overview:
##
##    freq *  mean *
## B 0.652 1 0.866 2
## A 0.544 0 0.594 0
## (* = original tree)
plot(st)

If I were to do some P-hacking I would combine the two features in a single LMM which leads to a somewhat better log-likelihood:

m <- lmer(Logit ~ factor(A >= 39.5) + B + (1 | AREA), data = d)
logLik(tr)
## 'log Lik.' -42.38439 (df=7)
logLik(m)
## 'log Lik.' -40.29202 (df=5)

But the corresponding summary (see below) also shows the correlation between the two variables. Hence, data-driven determination of what is going on is hard on this moderately-sized sample. Small changes in the model specification and/or the data will prefer either A or B or both...

summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Logit ~ factor(A >= 39.5) + B + (1 | AREA)
##    Data: d
##
## REML criterion at convergence: 80.6
##
## Scaled residuals:
##      Min       1Q   Median       3Q      Max
## -2.11309 -0.42937  0.08232  0.55221  1.95544
##
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  AREA     (Intercept) 0.2522   0.5022
##  Residual             0.2025   0.4500
## Number of obs: 50, groups:  AREA, 13
##
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            -2.3493     0.6021  -3.902
## factor(A >= 39.5)TRUE  -0.8524     0.3918  -2.175
## B                       3.0293     1.2582   2.408
##
## Correlation of Fixed Effects:
##             (Intr) f(A>=3
## f(A>=39.5)T -0.708
## B           -0.943  0.571