0

Following the model-based recursive partitioning in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6015941/ I want to replicate the following code:

sim_data <- function(n=2000){
  x1 <- rnorm(n)
x2 <- rbinom(n,1,0.3)
x3 <- runif(n)
x4 <- rnorm(n)
t <- rbinom(n,1,0.5)
z <- 1-x2+x1+2*(x1>=0)*x2*t-2*(x1<0)*x2*t
pr <- 1/(1+exp(-z))
y <- as.factor(rbinom(n,1,pr))
data.frame(x1,x3,x2=as.factor(x2),x4, t=factor(t,labels=c("C","A")),y,z)
}
dt <- sim_data()

dt.num = as.data.frame(sapply(dt, as.numeric))
dt.num$y <- dt.num$y-1    #only to convert outcome 1,2 into 0,1
mbase <- glm(y~t, data=dt.num,
             family = binomial())
round(summary(mbase)$coefficients,3)
library("model4you")
pmtr <- pmtree(mbase, zformula = ~. ,
               data = dt.num,
               control = ctree_control(minbucket = 250))

plot(pmtr, terminal_panel = node_pmterminal(pmtr,
                                            plotfun = binomial_glm_plot,
                                            confint = TRUE))

However, the following inexplicable error occurs:

Error in .Call.graphics(C_palette2, .Call(C_palette2, NULL)) : 
  invalid graphics state

I was looking for a solution to this problem in the post Persistent invalid graphics state error when using ggplot2. But the problem persists.

Any clue?

Thank you in advance

Achim Zeileis
  • 15,710
  • 1
  • 39
  • 49
vog
  • 770
  • 5
  • 11

1 Answers1

0

When I tried to replicate this, I got a different error:

plot(pmtr, terminal_panel = node_pmterminal(pmtr, plotfun = binomial_glm_plot, confint = TRUE))
## Waiting for profiling to be done...
## Error in plotfun(mod = list(coefficients = c(`(Intercept)` = -0.16839363929017,  : 
##   Plotting currently only works for models with a single factor covariate. 
##         We recommend using partykit or ggparty plotting functionalities!

The reason for this is that the panel function expects both the response and the treatment to be binary factors (as in dt). When you use binary numeric variables instead (as in dt.num) the model estimation in glm() leads to equivalent output but the plot() functionality is confused.

When I refit both the glm() and the pmtree() with dt rather than dt.num everything works as intended for me, yielding the following graphic:

pmtree plot

Achim Zeileis
  • 15,710
  • 1
  • 39
  • 49