11
$\begingroup$

Data precisions :

  • quotation is a dummy variable
  • minute count all the minutes within a day
  • temp is the temperature

Here is my code :

ctree <- ctree(quotation ~ minute + temp, data = visitquot)
print(ctree)

Fitted party:
[1] root
|   [2] minute <= 600
|   |   [3] minute <= 227
|   |   |   [4] temp <= -0.4259
|   |   |   |   [5] temp <= -2.3174: 0.015 (n = 6254, err = 89.7)
|   |   |   |   [6] temp > -2.3174
|   |   |   |   |   [7] minute <= 68: 0.028 (n = 4562, err = 126.3)
|   |   |   |   |   [8] minute > 68: 0.046 (n = 7100, err = 312.8)
|   |   |   [9] temp > -0.4259
|   |   |   |   [10] temp <= 6.0726: 0.015 (n = 56413, err = 860.5)
|   |   |   |   [11] temp > 6.0726: 0.019 (n = 39779, err = 758.9)
|   |   [12] minute > 227
|   |   |   [13] minute <= 501
|   |   |   |   [14] minute <= 291: 0.013 (n = 30671, err = 388.0)
|   |   |   |   [15] minute > 291: 0.009 (n = 559646, err = 5009.3)
|   |   |   [16] minute > 501
|   |   |   |   [17] temp <= 5.2105
|   |   |   |   |   [18] temp <= -1.8393: 0.009 (n = 66326, err = 617.1)
|   |   |   |   |   [19] temp > -1.8393: 0.012 (n = 355986, err = 4289.0)
|   |   |   |   [20] temp > 5.2105
|   |   |   |   |   [21] temp <= 13.6927: 0.014 (n = 287909, err = 3900.7)
|   |   |   |   |   [22] temp > 13.6927
|   |   |   |   |   |   [23] temp <= 14: 0.035 (n = 2769, err = 92.7)
|   |   |   |   |   |   [24] temp > 14: 0.007 (n = 2161, err = 15.9)
|   [25] minute > 600
|   |   [26] temp <= 1.6418
|   |   |   [27] temp <= -2.3366: 0.012 (n = 110810, err = 1268.1)
|   |   |   [28] temp > -2.3366: 0.014 (n = 584457, err = 7973.2)
|   |   [29] temp > 1.6418: 0.016 (n = 3753208, err = 57864.3)

Then I ploted the tree :

plot(ctree, type = "simple")

And here is a part of the output :

enter image description here

My questions are :

  1. In the first output from print(ctree), lets take the last line [29] temp > 1.6418: 0.016 (n = 3753208, err = 57864.3). What does the value 0.016 means ? is that a p-value ? And what does err = 57864.3 means ? It can't be a count of attribution error because it is a float number.
  2. I could not find anywhere a similar output that I have in the grey square. If someone knows how to interpret it. And how can a p-value be negative ?

Here same code but different output

Improve this question
$\endgroup$
5
  • $\begingroup$ I would suggest to start by trying the same process using party package instead of partykit. I think the type='simple' plot works better that way. What is your dummy variables' values? Is it binary, categorical? Is this a classification tree or regression tree? Would be good to see a summary of your 3 variables. I have a feeling that your dummy (output) variable is numeric but the model treats it as scale variable and not categorical. $\endgroup$
    – AntoniosK
    Commented Sep 6, 2015 at 12:36
  • $\begingroup$ Then I think that 0.016 you get is the average of the 0s and 1s. You can see it as the % of 1s (success classes), but the right way to do it is to use as.character instead of using factors. You have to let the model know that those 0s and 1s are labels and not real numbers. $\endgroup$
    – AntoniosK
    Commented Sep 6, 2015 at 12:56
  • $\begingroup$ The variable isnumeric and quotation is composed of 0 and 1. How can I know if it is a regression or categorical tree ? I will test right away with quotation as a factor and then with the party package. However, what I understand from your answer is that the current output isn't normal right ? $\endgroup$
    – Yohan Obadia
    Commented Sep 6, 2015 at 13:01
  • $\begingroup$ If your output variable is a scale variable the method recognises it and builds a regression tree. If your output is categorical the method will build a classification tree. There's also some hint based on your err values. If it's a classification tree those will be a missclasification %. $\endgroup$
    – AntoniosK
    Commented Sep 6, 2015 at 13:09
  • $\begingroup$ As for the plot(..., type = "simple") problem. I still need to check why this currently does not work as desired in partykit but will try to fix this soon. In the meantime, just do plot(as.simpleparty(ctree)) which generates the desired plot. (This is better than going back to the old party implementation...) $\endgroup$
    – Achim Zeileis
    Commented Sep 6, 2015 at 19:50

3 Answers 3

Reset to default
13
$\begingroup$

So the 0.016 that you are seeing is the average of quotation while err is just the SSE.

I'm old school and don't know how to exactly show this using the new partykit package (maybe @Achim could illustrate), but I can show you how this done using the older party package.

So let's create some reproducible example

library(partykit)
airq <- subset(airquality, !is.na(Ozone))
airct <- ctree(Ozone ~ ., data = airq)
print(airct)
# Model formula:
#   Ozone ~ Solar.R + Wind + Temp + Month + Day
# 
# Fitted party:
#   [1] root
# |   [2] Temp <= 82
# |   |   [3] Wind <= 6.9: 55.600 (n = 10, err = 21946.4)
# |   |   [4] Wind > 6.9
# |   |   |   [5] Temp <= 77: 18.479 (n = 48, err = 3956.0)
# |   |   |   [6] Temp > 77: 31.143 (n = 21, err = 4620.6)
# |   [7] Temp > 82
# |   |   [8] Wind <= 10.3: 81.633 (n = 30, err = 15119.0)
# |   |   [9] Wind > 10.3: 48.714 (n = 7, err = 1183.4)
# 
# Number of inner nodes:    4
# Number of terminal nodes: 5

Now let's dtach the partykit package and fit the same tree using party and calculate the same values

detach("package:partykit", unload=TRUE)
library(party)
airct <- party::ctree(Ozone ~ ., data = airq)

t(sapply(unique(where(airct)), function(x) {
  n <- nodes(airct, x)[[1]]
  Ozone <- airq[as.logical(n$weights), "Ozone"]
  cbind.data.frame("Node" = as.integer(x), 
                   "n" = length(Ozone), 
                   "Avg."= mean(Ozone), 
                   "Variance"= var(Ozone), 
                   "SSE" = sum((Ozone - mean(Ozone))^2))  
}))

#      Node n  Avg.     Variance SSE     
# [1,] 5    48 18.47917 84.16977 3955.979
# [2,] 3    10 55.6     2438.489 21946.4 
# [3,] 6    21 31.14286 231.0286 4620.571
# [4,] 9    7  48.71429 197.2381 1183.429
# [5,] 8    30 81.63333 521.3437 15118.97

I think it is easy to see from the output which is which

What the code above does is basically extract the information of Ozone information out of the fitted tree per each inner node and calculate the relevant statistics

As per the P-values, this seems like some type of printing bug, here's an example how you can calculate the P-values in the inner nodes

terNodes <- unique(where(airct))
setdiff(1:max(terNodes), terNodes)

sapply(setdiff(1:max(terNodes), terNodes), function(x) {
          n <- nodes(airct, x)[[1]]
          pvalue <- 1 - nodes(airct, x)[[1]]$criterion$maxcriterion
          plab <- ifelse(pvalue < 10^(-3),
                         paste("p <", 10^(-3)),
                         paste("p =", round(pvalue, digits = 3)))
          c("Node" = x, "P-value" = plab)
})

#         [,1]        [,2]        [,3]        [,4]       
# Node    "1"         "2"         "4"         "7"        
# P-value "p < 0.001" "p = 0.002" "p = 0.003" "p = 0.003"

As a side note, if quotation is a dummy variable, you should probably fit a classification tree (i.e. convert it to a factor) rather a regression tree like you are doing

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ I can't upvote your answer yet due to lack of reputation but thanks. Your answer combined with @Antoniosk helped me greatly ! $\endgroup$
    – Yohan Obadia
    Commented Sep 6, 2015 at 14:36
  • $\begingroup$ Regarding that negative P-value you are seing. It looks like a bug to me as P-value can't be negative by definition. I'm sure @Achim will have something to add on that matter. $\endgroup$
    – David Arenburg
    Commented Sep 6, 2015 at 14:56
  • $\begingroup$ Anyway, also added an illustration regarding how to calculate p-values in the inner nodes- see edit. $\endgroup$
    – David Arenburg
    Commented Sep 6, 2015 at 15:11
  • $\begingroup$ To obtain the node-specific summaries, the simplest solution (IMHO) is to just do a tapply() on the response by the corresponding terminal nodes, computing any function you desire, e.g., tapply(airq$Ozone, predict(airct, type = "node"), function(y) c("n" = length(y), "Avg." = mean(y), "Variance" = var(y), "SSE" = sum((y - mean(y))^2))). You can also do a subsequent do.call("rbind", ...) or use some other kind of aggregation... $\endgroup$
    – Achim Zeileis
    Commented Sep 6, 2015 at 19:58
  • $\begingroup$ To extract the p-values, you can easily extract these in the new partykit version. To obtain the p-values from all tests carried out, just do library("strucchange") and then sctest(airct). From this you can easily get the minimum or any other summary you desire. Furthermore, youcan also just extract the minimal p-value stored in each node (if any), by nodeapply(airct, ids = nodeids(airct), FUN = function(n) info_node(n)$p.value). $\endgroup$
    – Achim Zeileis
    Commented Sep 6, 2015 at 20:02
6
$\begingroup$

As suggested by @DavidArenburg, I'm collecting my comments here in another answer for easier references.

(1) If your response is binary, you need to turn it into a factor. Otherwise the inference employed when growing the tree is not what it should be, and also the predictions, visualizations, and error measures are not those intended for classification. See the replies by @DavidArenburg and @AntoniosK for further examples. In general: Also the explanatory variables need to have appropriate classes (numeric, factor, ordered factor) to be processed correctly when growing the tree.

(2) The plot(..., type = "simple") currently does not work as desired - in other words this is a bug. We will fix the partykit package in due course. For the moment you can easily work around it by using plot(as.simpleparty(...)). As a reproducible example:

library("partykit")
data("PimaIndiansDiabetes", package = "mlbench")
ct <- ctree(diabetes ~ ., data = PimaIndiansDiabetes)
plot(as.simpleparty(ct))

simpleparty

(3) The $criterion table erroneously shown currently in the type = "simple" plot contains the test statistics and corresponding log 1-p values from the conditional inference conducted in each node. The log rather than the p-value is used because it is numerically much more stable when used for comparisons, computing the minimal value, etc. Note that the p-values can become extremely small when significant. As an example consider the test results corresponding to the first node from the tree above:

nodeapply(ct, ids = 1, function(n) info_node(n)$criterion)
## $`1`
##                pregnant       glucose   pressure   triceps      insulin
## statistic  3.631413e+00  5.117841e+00  1.1778530  1.455334  2.570457503
## p.value   -6.380290e-09 -2.710725e-37 -0.5937987 -0.313498 -0.002398554
##                    mass      pedigree           age
## statistic  4.185236e+00  3.143294e+00  3.774507e+00
## p.value   -4.181295e-15 -1.180135e-05 -3.262459e-10

(4) To extract the actual p-values (without log), there is an extractor function sctest() (for structural change test) which is used by the strucchange package and also for the parameter instability tests in the mob() trees. In the example above:

library("strucchange")
sctest(ct, node = 1)
##               pregnant  glucose  pressure   triceps     insulin         mass
## statistic 3.776615e+01 166.9745 3.2473947 4.2859164 13.07180346 6.570902e+01
## p.value   6.380290e-09   0.0000 0.4477744 0.2691142  0.00239568 4.218847e-15
##               pedigree          age
## statistic 2.318009e+01 4.357601e+01
## p.value   1.180128e-05 3.262459e-10

Note that one p-value (glucose) becomes zero while its log(1 - p) is very close to zero but not quite. To only extract the minimal p-values from the selected partitioning variable (if any), you can again use the nodeapply() function to obtain it from each node's $info:

nodeapply(ct, ids = nodeids(ct), function(n) info_node(n)$p.value)
## $`1`
## glucose 
##       0 
## 
## $`2`
##          age 
## 6.048661e-07 
## 
## $`3`
##        mass 
## 0.001169778 
## 
## ...

(5) If you want to summarize the response in each terminal node, the easiest solution (IMHO) is to do a simple tapply() on the response variable by the node groups, supplying any summary function you want. For a simpler overview, you can also rbind() this etc. For example:

tab <- tapply(PimaIndiansDiabetes$diabetes, predict(ct, type = "node"),
  function(y) c("n" = length(y), 100 * prop.table(table(y))))
do.call("rbind", tab)
##      n      neg        pos
## 5  144 99.30556  0.6944444
## 6    7 85.71429 14.2857143
## 7  120 82.50000 17.5000000
## 8  214 66.82243 33.1775701
## 11  53 79.24528 20.7547170
## 12 108 39.81481 60.1851852
## 13 122 19.67213 80.3278689
Improve this answer
$\endgroup$
3
$\begingroup$

As an addition to @DavidArenburg 's great answer I'll show you the difference in outputs between regression and classification trees

library(partykit)

# data
airquality = data.frame(airquality)

# create a numeric binary variable as dependent variable
airquality$OzoneClass = 0
airquality$OzoneClass[airquality$Ozone>=34] =1

# regression tree with scale dependent variable
airq <- subset(airquality, !is.na(Ozone))
airct <- ctree(OzoneClass ~ Temp, data = airq)
print(airct)

# Model formula:
#   OzoneClass ~ Temp
# 
# Fitted party:
#   [1] root
# |   [2] Temp <= 82
# |   |   [3] Temp <= 77: 0.096 (n = 52, err = 4.5)
# |   |   [4] Temp > 77: 0.519 (n = 27, err = 6.7)
# |   [5] Temp > 82: 0.973 (n = 37, err = 1.0)
# 
# Number of inner nodes:    2
# Number of terminal nodes: 3




# create a categorical binary variable as dependent variable
airquality$OzoneClass = 0
airquality$OzoneClass[airquality$Ozone>=34] =1
airquality$OzoneClass = as.factor(airquality$OzoneClass)

# classification tree
airq <- subset(airquality, !is.na(Ozone))
airct <- ctree(OzoneClass ~ Temp, data = airq)
print(airct)

# Model formula:
#   OzoneClass ~ Temp
# 
# Fitted party:
#   [1] root
# |   [2] Temp <= 82
# |   |   [3] Temp <= 77: 0 (n = 52, err = 9.6%)
# |   |   [4] Temp > 77: 1 (n = 27, err = 48.1%)
# |   [5] Temp > 82: 1 (n = 37, err = 2.7%)
# 
# Number of inner nodes:    2
# Number of terminal nodes: 3

Note how (a) the values of the trees change from an average value of 0s and 1s (regression tree) to the most likely class 0/1 (classification tree) and (b) the err values from a number (error) become a percentage (missclassification).

Improve this answer
$\endgroup$
0

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.