Details of test statistic used for conditional inference regression trees

Question

Conditional inference regression trees implemented in the partykit R-library use the following statistic for determining if a split will be carried out or not:

$$ T_j(L_n, w) = vec(\sum w_i g_j(X_{ji}) h(Y_i, (Y_1,...,Y_n)^T))$$

(See Variable selection and stopping criteria in ctree.pdf)

According to https://stats.stackexchange.com/a/342356/161138 this independence test can be replicated with the coin library which works for me. It can be further broken down into its parts as follows:

library(coin)
it <- independence_test(Petal.Length ~ Sepal.Length,
  data = iris, teststat = "quadratic")

(statistic(it, "linear") - expectation(it)) %*% 
   solve(covariance(it)) %*% 
   t(statistic(it, "linear") - expectation(it))

I have simplified the example to just one independent variable Petal.Length ~ Sepal.Length, there are two in the original question Petal.Length + Petal.Width ~ Sepal.Length.

Now, I am trying to reproduce all the steps in my simplified case:

• statistic(it, "linear") is identical to sum(iris$Petal.Length * iris$Sepal.Length)
• expectation(it) is identical to sum(iris$Petal.Length) * sum(iris$Sepal.Length) / length(iris$Petal.Length)

What is

• covariance(it) in terms of iris$Petal.Length & iris$Sepal.Length?

Ultimately I want to write the first equation for my specific case of just one independent variable. I believe it comes down to $$ \frac{(t-\mu)^2}{\sigma^2} $$ where

• $ t = \sum (a \times b) $,
• $\mu = \frac{\sum a \times \sum b}{n}$ &
• $\sigma^2 = ?$.

I am looking for a formula for $\sigma$.

I guess this is described here: https://cran.r-project.org/web/packages/coin/vignettes/MAXtest.pdf, but I can't figure out how to rewrite $\sigma$ for my very simple case.

I thought this is the answer, but it is wrong

(n/(n-1)) *sum(iris$Petal.Length * iris$Sepal.Length) - (1/(n-1)) * sum(iris$Petal.Length) * sum(iris$Sepal.Length)

which would mean $\sigma^2$ is (but it is not)

$$ \frac{n}{n-1} \sum (a \times b) - \frac{1}{n-1} \times \sum a \times \sum b $$

Answer 1

The equations are easier to follow if you use the same notation. Also, the equations become simpler because both transformations are the identity, i.e., $h(Y_i) = Y_i$ and $g(X_i) = X_i$. Hence, the transposes don't matter and all the matrix multiplications and Kronecker products are just scalar products as well.

So the two variables are the following and we don't need any transformations here:

X <- iris$Sepal.Length
Y <- iris$Petal.Length

Then we need the mean and variance of Y:

E <- sum(Y)/N
V <- sum((Y - E)^2)/N

Then the mean and variance of the linear statistic are:

mu <- E * sum(X)
Sigma <- N/(N - 1) * V * sum(X^2) - V/(N - 1) * sum(X)^2

And then we can compute the linear statistic as well as the quadratic statistic:

T <- sum(X * Y)
(T - mu) * Sigma^(-1) * (T - mu)
## [1] 113.2332

Compare that with the statistic in it:

statistic(it)
## [1] 113.2332

And the building blocks are also equal:

all.equal(T, c(statistic(it, "linear")))
## [1] TRUE
all.equal(mu, c(expectation(it)))
## [1] TRUE
all.equal(Sigma, c(covariance(it)))
## [1] TRUE

With scalar transformations virtually nothing changes except that in the beginning X and Y get transformed.

The computations become more involved when we use multivariate transformation because then all the matrix operations matter. However, for the most important case of multivariate transformation - namely dummy variables for all factor levels - the computations become relatively easy again...