#r #decision-tree #party
#r #дерево решений #Вечеринка
Вопрос:
Я использую mob() функцию из partykit пакета, и у меня возникают некоторые проблемы при анализе полученной модели.
Моя главная цель — приблизительно проверить, насколько большим должен быть размер выборки, чтобы определить реальную структуру процесса генерации данных (DGP) при наличии перерывов.
Приведенный ниже код выполняет моделирование данных Montecarlo с разрывами и пытается определить, был ли разрыв зафиксирован тестом M-fluctuation или нет.
Более конкретно, я хочу подсчитать, сколько раз по сравнению с общим количеством симуляций ( nreps ) модель фактически фиксирует DGP, при условии фиксированного размера выборки ( N ), чтобы получить представление о том, сколько данных мне нужно для захвата реального DGP.
В конце кода вы можете увидеть список, который я получаю из своих симуляций. Проблема в том, что я не могу восстановить информацию, отображаемую на консоли.
Кроме того, у меня есть некоторые сомнения относительно того, как сделать подсчет «правильно идентифицированных моделей». На данный момент я имею в виду считать положительным, если модель имеет разрыв в правильную переменную ( z1 ) в указанной области ( z1>0 ) с некоторым допуском в область разрыва, например, если разрыв равен z1>0.1 или z1>-0.4 он также действителен как положительный для меня. Поэтому существует ли какой-либо простой способ подсчета моделей, обладающих указанными характеристиками?
Я надеюсь, что мой пример достаточно понятен, чтобы вы могли мне помочь. Заранее большое вам спасибо.
Имитировать модель
- Ингредиенты для создания DGP.
library("partykit")
library(data.table) ## optional, but what I'll use to coerce the list into a DT
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
#sample size
N <- 300
#coeficients
betas <- list()
betas$b0 <- 1
betas$b1_up <- 2.4
betas$b1_down <- 2
betas$b2_up <- 2.4
betas$b2_down <- 2
#mob() ingredients
ols_formula <- y ~ x1 x2 | z1 z2
# the ""0 "" ---> supress the 'double' interecept
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 x)}
- Функция, которая генерирует данные и подгоняет OLS с использованием алгоритма mob.
reg_simulation_mob <- function(...){
#data
data <- data.frame(
x1 = rnorm(N),
x2 = rnorm(N),
z1 = rnorm(N),
z2 = rnorm(N),
e = rnorm(N))
#dependent variable
data$y <- betas$b0 with(data, ifelse(z1>0,
betas$b1_up * x1 betas$b2_up * x2 ,
betas$b1_down * x1 betas$b2_down * x2 )
e )
#Estimate mob()-OLS
ols_mob <- mob(ols_formula,
data = data,
fit = ols)
# return(ols$coefficients)
return(ols_mob)
}
- Моделирование Montecarlo (только 2 испытания) описанной выше настройки.
# N repetitions
nreps <- 2
## Parallel version
results <- future_lapply(1:nreps, reg_simulation_mob, future.seed = 1234L)
Полученный результат
Как вы можете видеть ниже, в первом испытании ( results[[1]] ) модель находит правильный разрыв, но во втором ( results[[2]] ) он не может его найти.
> results
[[1]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 x2 | z1 z2
Fitted party:
[1] root
| [2] z1 <= 0.00029: n = 140
| x(Intercept) xx1 xx2
| 0.9597894 1.7552122 2.1360788
| [3] z1 > 0.00029: n = 160
| x(Intercept) xx1 xx2
| 0.9371795 2.4745728 2.5087608
Number of inner nodes: 1
Number of terminal nodes: 2
Number of parameters per node: 3
Objective function: 422.2329
[[2]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 x2 | z1 z2
Fitted party:
[1] root: n = 300
x(Intercept) xx1 xx2
1.015224 2.175625 2.200746
Number of inner nodes: 0
Number of terminal nodes: 1
Number of parameters per node: 3
Objective function: 422.3085
На рисунке ниже вы можете наблюдать структуру списка results , где я не могу найти информацию, отображаемую на консоли (т.Е. Количество узлов, параметры, пороговые значения и т. Д.)
Ответ №1:
Во-первых, я бы рекомендовал использовать lmtree() функцию, а не ваниль mob() . Первый работает быстрее, имеет лучшие настройки по умолчанию для печати и построения графиков и имеет больше возможностей для прогнозирования.
Во-вторых, я рекомендую вам обратиться vignette("partykit", package = "partykit") к, в котором объясняется, как party создаются объекты и какие классы и методы задействованы.
В-третьих, чтобы определить, какая переменная (если таковая имеется) использовалась для разделения в корневом узле, вероятно, представляет интерес извлечь результаты всех тестов нестабильности параметров. sctest() Для получения этого существует специальный метод (проверка структурных изменений):
library("strucchange")
sctest(results[[1]], node = 1)
## z1 z2
## statistic 4.072483e 01 6.1762164
## p.value 5.953672e-07 0.9153013
sctest(results[[2]], node = 1)
## z1 z2
## statistic 12.2810548 10.1944484
## p.value 0.2165527 0.4179998
Соответствующий partysplit объект для $split (если таковой имеется) в корне $node , вероятно, проще всего извлечь вручную:
results[[1]]$node$split
## $varid
## [1] 4
##
## $breaks
## [1] 0.0002853492
##
## $index
## NULL
##
## $right
## [1] TRUE
##
## $prob
## NULL
##
## $info
## NULL
##
## attr(,"class")
## [1] "partysplit"
results[[2]]$node$split
## NULL
Идентификатор переменной соответствует порядку переменных в:
names(results[[1]]$data)
## [1] "y" "x1" "x2" "z1" "z2"
Наконец, что касается вопроса, что оценивать: все зависит от определения правильной переменной для разделения. Если это сделано правильно, то оценки точек разделения быстро сходятся к истинным значениям, а затем оценки параметров также сходятся. См., Например, нашу недавнюю статью arXiv (https://arxiv.org/abs/1906.10179 ), который содержит более масштабное исследование моделирования, а также предоставляет код репликации.
Поэтому, как правило, я либо оцениваю вероятность выбора правильной переменной при первом разделении. Или, в качестве альтернативы, я смотрю на RMSE оцененных и истинных коэффициентов для каждого наблюдения.
Обновление: за пределами корневого узла вы можете использовать nodeapply() для извлечения информации из различных узлов. Однако я не рекомендую оценивать все разбиения, потому что становится все труднее сопоставить, какое оцененное разделение соответствует какому из истинных разбиений. Вместо этого мы часто оцениваем, насколько подобранный раздел сравнивается с истинным разделом, например, используя скорректированный индекс Rand. Опять же, вы можете найти пример для в упомянутой выше статье arXiv.
Комментарии:
1. Большое спасибо за ваш любезный ответ, профессор @AchimZeileis. Я прочитал то, что вы мне предложили, и это очень помогло мне лучше понять
partyобъекты. Однако у меня есть последний вопрос, и я воспользовался свободой редактирования своего вопроса с помощью расширенного примера, который должен проиллюстрировать мою проблему и желаемый результат. Еще раз спасибо.2. Я добавил несколько обновленных комментариев в свой ответ. Однако обратите внимание, что такие бесконечные вопросы-спагетти могут быть полезны для вас, но (а) обычно на них трудно ответить и (б) они менее полезны для других читателей. Поэтому рекомендуется задавать целенаправленные вопросы о конкретной проблеме, на которые, вероятно, есть точный ответ. Более подробную информацию см. В руководстве по вопросам.
3. Спасибо за обновление в вашем ответе, профессор @AchimZeileis. Кроме того, приношу свои извинения за беспорядок в последнем редактировании моего первоначального вопроса. Я изменил его, чтобы ваш ответ оставался действительным, и, кроме того, я добавил расширенный ответ на свой собственный вопрос, который, я надеюсь, может быть полезен другим пользователям в будущем. Большая часть кода и рисунков были взяты (и процитированы) из упомянутой вами ссылки. Еще раз спасибо.
4. Хорошо, спасибо. Это дольше и сложнее, чем то, что toy обычно публикует в одном вопросе / ответе на SO. В будущем было бы лучше разбить такие обсуждения на несколько вопросов. Для этого вопроса я бы рекомендовал оставить все как есть, но принять один из двух ответов, чтобы вопрос был отмечен решенным здесь на SO.
Ответ №2:
Этот ответ основан на ссылке, которую профессор @AchimZeileis предоставил в своей статье (https://arxiv.org/abs/1906.10179 ), и он посвящен второй части моего первоначального вопроса, который касался вопроса: как подсчитать правильно указанные модели (деревья)?
Короткий ответ.
В статье проблема разделена на два разных типа процесса генерации данных (DGP). В первом случае данные имеют только один разрыв в одной переменной (случай «тупика«), и авторы подсчитывают количество правильно идентифицированных моделей на основе количества раз, когда тест M-fluctuation правильно идентифицировал переменную, которая создавала разрыв (только один реальный разрыв в одной переменнойи 9 шумных кандидатов для разделения переменной без разрыва). Второй DGP был моделью с двумя разрывами (случай «дерева»), и они использовали скорректированный индекс Ранда (ARI) для оценки производительности модели как показателя сходства реального дерева с предсказанным.
Очень длинный ответ.
Давайте разберем ARI на 6 различных иллюстративных возможных деревьев, которые могут быть получены при разных размерах выборки. Используемый здесь код в значительной степени основан на дополнительном материале статьи, рекомендованной @AchimZeileis.
Процесс генерации данных: древовидная структура
Реальный dgp имеет 2 перерыва, как показано на рисунке ниже. Первая генерируется переменной z2 , а вторая z1 — . Во фрагменте кода ниже дельта равна 1. Пороговое значение для первого разрыва (в зависимости от z2 ) равно 0,3, а пороговое значение для второго разрыва (в зависимости от z1 ) равно -0,3 (значения можно увидеть в объекте xi = c(-0.3, 0.3) )
#function from https://arxiv.org/src/1906.10179v1/anc
dgp_tree <- function(nobs = 1000, delta = 1, xi = c(-0.3, 0.3),
sigma = 1, seed = 7,
changetype = "abrupt",
variation = "all",
beta0 = NULL, beta1 = NULL)
{
# check input values
if(variation != "all") stop("variation can only be set to 'all' in dgp_tree")
if(changetype != "abrupt") stop("changetype can only be abrupt in dgp_tree")
if(!is.null(beta0) | !is.null(beta1)) warning("values for beta0 or beta1 are ignored since variation='all' for dgp_tree")
set.seed(seed)
if(length(xi)==1){
xi1 <- xi2 <- xi
} else {
xi1 <- xi[1]
xi2 <- xi[2]
}
z1 <- runif(nobs,-1,1)
z2 <- runif(nobs,-1,1)
z3 <- rnorm(nobs, 0, 1)
z4 <- rnorm(nobs, 0, 1)
z5 <- rnorm(nobs, 0, 1)
z6 <- rnorm(nobs, 0, 1)
z7 <- rnorm(nobs, 0, 1)
z8 <- runif(nobs, -1, 1)
z9 <- runif(nobs, -1, 1)
z10 <- runif(nobs, -1, 1)
id <- numeric(length(z1))
x <- runif(nobs, min = -1, max = 1)
beta0 <- delta * (-1)^(z1<xi1) * 0^(z2<xi2)
beta1 <- delta * (-1)^(z2>=xi2)
id <- 1 (z2>=xi2) (z2>=xi2)*(z1>=xi1)
mu <- beta0 beta1 * x
y <- rnorm(nobs, mu, sigma)
d <- data.frame(y = y, x = x,
z1 = z1, z2 = z2, z3 = z3, z4 = z4, z5 = z5, z6 = z6, z7 = z7, z8 = z8, z9 = z9, z10 = z10,
beta0 = beta0, beta1 = beta1, mu = mu, sigma = rep.int(sigma, times = length(y)), id = id)
return(d)
}
Скорректированный индекс Rand
Среди функций, включенных в статью, есть одна для вычисления ARI, и она указана ниже для использования в следующих примерах. Это почти точно напоминает букву за буквой обозначение, используемое здесь.
# function to compute adjusted Rand Index from https://arxiv.org/src/1906.10179v1/anc
adj_rand_index <- function(x, y) {
tab <- table(x, y)
a <- rowSums(tab)
b <- colSums(tab)
M <- sum(choose(tab, 2))
N <- choose(length(x), 2)
A <- sum(choose(a, 2))
B <- sum(choose(b, 2))
c(ARI = (M - (A * B) / N) / (0.5 * (A B) - (A * B) / N))
}
Моделирование данных
library(partykit)
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
ols_formula <- y ~ x | z1 z2 z3 z4 z5 z6 z7 z8 z9 z10
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 x)}
sim_ari <- function(n){
tree_data <- dgp_tree(nobs = n)
ols_mob <- mob(ols_formula,
data = tree_data,
fit = ols)
prednode <- predict(ols_mob ,
type = "node")
cross_table <- table(prednode,tree_data$id)
ari <- adj_rand_index(prednode,
tree_data$id)
print(n)
print(ari)
return(
list(
ols_mob = ols_mob,
cross_table = cross_table,
ari=ari,
data = tree_data)
)
}
n_levels <- c(55, ## no break
87, ## only one break
123, ## Correct structure, but poor performance
199, ## Nested break in second leaf
667, ## Additional break in first leaf
5000 ## Perfect model
)
ari <- future_lapply(n_levels, sim_ari, future.seed = 1234L)
A Tale of Six Trees.
The following six cases are analyzed in terms of how the ARI can accurately capture the degree of similarity between the correct and the estimated tree. The key to compare the trees is id , which shows which leaf each observation should belong in the tree. For example, if an observation has an id value of 1 , it meets the requirements assigned to node number 2 in the figure above. On the other hand, if id is equal to 2, the observation should be assigned to node number 4 in the same picture. Finally, if id is equal to 3, it is assigned to node number 5. You can check this reasoning in the following line id <- 1 (z2>=xi2) (z2>=xi2)*(z1>=xi1)
The First Tree (n=55): No breaks
The first case analyzed corresponds when no breaks are identified. Here, in this case, the ARI is equal to 0.
##### First Tree (n=55): No break ####
ari[[1]][[1]]
## Fitted party:
## [1] root: n = 55
## x(Intercept) xx
## -0.01309586 0.39291089
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
## Number of parameters per node: 2
## Objective function: 95.58631
Здесь интересно отметить, что все наблюдения назначаются корневому узлу. Поэтому при перекрестной настройке прогнозируемых узлов prednode_1 мы видим, что все возможные id значения принадлежат корневому узлу [1] прогнозируемого дерева (в основном, потому что другого варианта нет.) Используя функцию adj_rand_index() , вы можете проверить, что это приводит к ARI, равному 0.
#data first tree (n=55)
data_1 <- ari[[1]][[4]]
#predicted node first iteration
data_1$prednode_1 <- predict(ari[[1]][[1]], type = "node")
#Cross table
with(data_1, table(prednode_1 ,id))
## id
## prednode_1 1 2 3
## 1 37 7 11
#adj_rand_index
ari[[1]][[3]]
Второе дерево (n = 87): выявлен только один разрыв
Этот случай интересен тем, что он частично идентифицирует структуру дерева (т. Е. Разрыв z1 отсутствует).
##### Second Tree (n=87): Extra partition in node[5] ####
ari[[2]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29288: n = 57
# | x(Intercept) xx
# | 0.133293 1.082701
# | [3] z2 > 0.29288: n = 30
# | x(Intercept) xx
# | 0.2598309 -1.8014133
#
# Number of inner nodes: 1
# Number of terminal nodes: 2
# Number of parameters per node: 2
# Objective function: 122.0116
Кроме того, мы можем проверить, что при перекрестной настройке прогнозируемых и реальных узлов мы видим, что некоторые наблюдения соответствуют критериям даже в этом несовершенном дереве. Это означает, что есть 57 наблюдения, правильно назначенные первому узлу и 9 которые были правильно назначены второй ветви. Наконец, 30 где неправильно назначено, потому что последний узел вообще не был идентифицирован. Это приводит к ARI, равному 0.8577366 , что является огромным улучшением по сравнению с первым деревом.
#data second iteration (n=87)
data_2 <- ari[[2]][[4]]
#predicted node first iteration
data_2$prednode_2 <- predict(ari[[2]][[1]], type = "node")
#Cross table
with(data_2, table(prednode_2 ,id))
# id
# prednode_2 1 2 3
# 2 57 0 0
# 3 1 9 20
#adj_rand_index
ari[[2]][[3]]
# > ari[[2]][[3]]
# ARI
# 0.8577366
Третье дерево (n = 123): правильная структура, но низкая производительность
Этот случай интересен тем, что он восстанавливает реальную древовидную структуру, но имеет худшую производительность, чем последние три, которые лишь частично идентифицировали его структуру.
##### Third Tree (n=123): Correct structure but poor performance ####
ari[[3]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.07319: n = 60
# | x(Intercept) xx
# | -0.1723388 1.1071878
# | [3] z2 > 0.07319
# | | [4] z1 <= -0.35485: n = 22
# | | x(Intercept) xx
# | | -0.7166565 -0.6791717
# | | [5] z1 > -0.35485: n = 41
# | | x(Intercept) xx
# | | 0.7096033 -0.8605967
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 156.4397
Ниже мы можем видеть, что когда мы сопоставляем прогнозируемые и реальные узлы, мы видим, что 16 ( 10 6 ) наблюдения были неправильно классифицированы, и этот факт приводит к ARI of 0.6117612 .
#data third iteration (n=123)
data_3 <- ari[[3]][[4]]
#predicted node first iteration
data_3$prednode_3 <- predict(ari[[3]][[1]], type = "node")
#Cross table
with(data_3, table(prednode_3 ,id))
# id
# prednode_3 1 2 3
# 2 60 0 0
# 4 6 16 0
# 5 10 0 31
#adj_rand_index
ari[[3]][[3]]
# > ari[[3]][[3]]
# ARI
# 0.6117612
Четвертое дерево: Четвертое дерево (n = 199): дополнительный лист в узле [5]
Указанное здесь дерево отклонилось от оригинала, потому что у него есть дополнительный лист из node[5] , которого нет в реальных данных.
##### Forth Tree (n=199): Extra leaf at node[5] ####
ari[[4]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= -0.19806: n = 79
# | x(Intercept) xx
# | 0.06455217 1.51512672
# | [3] z2 > -0.19806
# | | [4] z1 <= -0.27127: n = 44
# | | x(Intercept) xx
# | | -0.4863122 -0.3860951
# | | [5] z1 > -0.27127
# | | | [6] z2 <= 0.17481: n = 23
# | | | x(Intercept) xx
# | | | -0.1335096 0.2046050
# | | | [7] z2 > 0.17481: n = 53
# | | | x(Intercept) xx
# | | | 1.0868488 -0.0290925
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 282.6727
Here the count of the crosstable of the real and predicted nodes is interesting because nodes [6] and [7] are inexistent in the real data, but they are getting observations that should, for example, be assigned to the node [1] ( 23 and 7 observations respectively.) This misallocation diminished the ARI index down to 0.4649789 .
#data forth iteration (n=199)
data_4 <- ari[[4]][[4]]
#predicted node first iteration
data_4$prednode_4 <- predict(ari[[4]][[1]], type = "node")
#Cross table
with(data_4, table(prednode_4 ,id))
# id
# prednode_4 1 2 3
# 2 79 0 0
# 4 16 27 1
# 6 23 0 0
# 7 7 0 46
#adj_rand_index
ari[[4]][[3]]
# ARI
# 0.4649789
The fifth Tree (n=667): Extra leaf at node [2]
This is another example of a tree with an incorrect structure where an extra leaf (based on a partition on z5 which is incorrect(!)) is attached to the node [2] .
##### Fifth Tree (n=667): Extra leaf at node[2] ####
ari[[5]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.28476
# | | [3] z5 <= 0.76285: n = 322
# | | x(Intercept) xx
# | | -0.1322881 0.9535337
# | | [4] z5 > 0.76285: n = 96
# | | x(Intercept) xx
# | | 0.1686863 1.3878776
# | [5] z2 > 0.28476
# | | [6] z1 <= -0.32001: n = 89
# | | x(Intercept) xx
# | | -0.9139858 -0.7957158
# | | [7] z1 > -0.32001: n = 160
# | | x(Intercept) xx
# | | 0.7661154 -0.8656553
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 927.9088
The crosstable from predicted and correct nodes shows us that most of the observations ( 322 ) that, in reality, belong to the first node [1] were assigned to the predicted node [3] . Finally, this poor structure leads to an ARI of 0.6932132 .`
#data third iteration (n=667)
data_5 <- ari[[5]][[4]]
#predicted node first iteration
data_5$prednode_5 <- predict(ari[[5]][[1]], type = "node")
#Cross table
with(data_5, table(prednode_5 ,id))
# id
# prednode_5 1 2 3
# 3 322 0 0
# 4 96 0 0
# 6 0 89 0
# 7 3 3 154
#adj_rand_index
ari[[5]][[3]]
# ARI
# 0.6932132
The sixth tree (n=5000): The Golden Tree!
This final tree recovers the data perfectly, both in the tree structure and in allocating the observation to each leaf.
##### Sixth Tree (n=5000): Extra leaf at node[2] ####
ari[[6]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29971: n = 3187
# | x(Intercept) xx
# | -0.008719923 1.022232280
# | [3] z2 > 0.29971
# | | [4] z1 <= -0.30286: n = 609
# | | x(Intercept) xx
# | | -0.9488846 -0.9813765
# | | [5] z1 > -0.30286: n = 1204
# | | x(Intercept) xx
# | | 1.0281410 -0.9565637
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 6992.848
Здесь мы можем видеть из перекрестной таблицы прогнозируемых и реальных узлов, что она идеально распределяет каждое наблюдение там, где оно принадлежит, что приводит к ARI, равному 1 .
#data sixt iteration (n=5000)
data_6 <- ari[[6]][[4]]
#predicted node first iteration
data_6$prednode_6 <- predict(ari[[6]][[1]], type = "node")
#Cross table
with(data_6, table(prednode_6 ,id))
# id
# prednode_6 1 2 3
# 2 3187 0 0
# 4 0 609 0
# 5 0 0 1204
#adj_rand_index
ari[[6]][[3]]
# ARI
# 1
Выводы.
Некоторые важные выводы можно извлечь из иллюстрации выше.
1.- ARI полезен для оценки степени сходства прогнозируемого дерева, которое может иметь структуру, сильно отличающуюся от реального дерева в процессе генерации данных.
2. — Восстановление правильной структуры дерева не приводит к ARI, равному единице.
3. — Неправильные деревья не обязательно будут ARI равными нулю.
Окончательное моделирование.
В заключение, вот небольшая симуляция, чтобы увидеть, как ведет себя индекс ARI при увеличении размера выборки.
### Final simulation
n_levels <-seq(from= 10, to= 2000, by= 5)
ari <- lapply(n_levels, sim_ari)
ari_models<- function(i){
ari <- ari_sim[[i]]$ari
n <- nobs(ari_sim[[i]]$ols_mob)
return(
list(ari = ari , n= n )
)
}
ari_n_list <- lapply(1:length(ari_sim), ari_models)
df <- data.frame(matrix(unlist(ari_n_list), nrow=length(ari_n_list), byrow=T))
colnames(df) <- c("ARI" , "N")
library(ggplot2)
ggplot(df, aes(N))
geom_line(aes(y = ARI, colour = "ARI"))
