Running 'backticks.R' [156s]
  Running 'cost.R' [1s]
  Comparing 'cost.Rout' to 'cost.Rout.save' ... OK
  Running 'cptest.R' [154s]
  Running 'minus_in_formula.R' [159s]
  Running 'priors.R' [166s]
  Running 'rescale.R' [1s]
  Comparing 'rescale.Rout' to 'rescale.Rout.save' ... OK
  Running 'testall.R' [1s]
  Comparing 'testall.Rout' to 'testall.Rout.save' ... OK
  Running 'treble.R' [176s]
  Running 'treble2.R' [171s]
  Running 'treble3.R' [176s]
  Running 'treble4.R' [0s]
  Comparing 'treble4.Rout' to 'treble4.Rout.save' ... OK
  Running 'usersplits.R' [162s]
  Running 'xpred1.R' [0s]
  Comparing 'xpred1.Rout' to 'xpred1.Rout.save' ... OK
  Running 'xpred2.R' [164s]
Running the tests in 'tests/backticks.R' failed.
Complete output:
  > ## allow backticks in rpart.matrix: see
  > ## https://stat.ethz.ch/pipermail/r-help/2012-May/314081.html
  > 
  > set.seed(10)
  > library(rpart)
  > Iris <- iris
  > names(Iris) <- sub(".", " ", names(iris), fixed=TRUE)
  > rpart(Species ~ `Sepal Length`, data = Iris)
  n= 150 
  
  node), split, n, loss, yval, (yprob)
        * denotes terminal node
  
  1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)  
    2) Sepal Length< 5.45 52   7 setosa (0.86538462 0.11538462 0.01923077) *
    3) Sepal Length>=5.45 98  49 virginica (0.05102041 0.44897959 0.50000000)  
      6) Sepal Length< 6.15 43  15 versicolor (0.11627907 0.65116279 0.23255814) *
      7) Sepal Length>=6.15 55  16 virginica (0.00000000 0.29090909 0.70909091) *
  > 
  > proc.time()
     user  system elapsed 
     0.21    0.01    0.21 
Running the tests in 'tests/cptest.R' failed.
Complete output:
  > library(rpart)
  > mystate <- data.frame(state.x77, region=factor(state.region))
  > names(mystate) <- c("population","income" , "illiteracy","life" ,
  +        "murder", "hs.grad", "frost",     "area",      "region")
  > #
  > # This little test came out of a query that cp did not scale
  > #  with the data size.  It does
  > #
  > # tdata = 20 copies of "mystate"  
  > # trees with tdata and trees with mystate should be the same (they are)
  > #  except for the n's
  > set.seed(10)
  > 
  > tdata <- rbind(mystate, mystate, mystate, mystate, mystate)
  > tdata <- rbind(tdata, tdata, tdata, tdata)
  > tfit1 <- rpart(income ~ population + illiteracy + murder + hs.grad + region,
  +                data = mystate, method = "anova", xval=0, cp=.089)
  > tfit2 <- rpart(income ~ population + illiteracy + murder + hs.grad + region,
  +                data = tdata, method='anova', xval=0, cp=.089, 
  +                minsplit=400, minbucket=140)
  > 
  > all.equal(tfit1$splits[,-1], tfit2$splits[,-1])
  [1] TRUE
  >  
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.06    0.21 
Running the tests in 'tests/minus_in_formula.R' failed.
Complete output:
  > ## https://github.com/bethatkinson/rpart/issues/7
  > 
  > library(rpart)
  > mtcars2 <- mtcars
  > 
  > mtcars2$gear <- factor(mtcars2$gear)
  > mtcars2$carb <- factor(mtcars2$carb)
  > 
  > set.seed(10)
  > rp1 <- rpart(mpg ~ . - gear, data = mtcars2)
  > set.seed(10)
  > rp2 <- rpart(mpg ~ ., data = mtcars2[names(mtcars2) != "gear"])
  > 
  > all.equal(rp1[setdiff(names(rp1), c("call", "terms"))],
  +           rp2[setdiff(names(rp2), c("call", "terms"))], check.attributes=FALSE)
  [1] TRUE
  > 
  > 
  > set.seed(10)
  > rp3 <- rpart(mpg ~ . - gear - carb, data = mtcars2)
  > set.seed(10)
  > rp4 <- rpart(mpg ~ ., data = mtcars2[setdiff(names(mtcars2), c("gear", "carb"))])
  > 
  > all.equal(rp3[setdiff(names(rp3), c("call", "terms"))],
  +           rp4[setdiff(names(rp4), c("call", "terms"))], check.attributes=FALSE)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.25    0.07    0.26 
Running the tests in 'tests/priors.R' failed.
Complete output:
  > #
  > # A hard-core test of losses and priors
  > #  Simple data set where I know what the answers must be
  > #
  > library(rpart)
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > dummy <- c(3,1,4,1,5,9,2,6,5,3,5,8,9,7,9)/5
  > pdata <- data.frame(y=factor(rep(1:3, 5)),
  +                     x1 = 1:15,
  +                     x2 = c(1:6, 1:6, 1:3),
  +                     x3 = (rep(1:3, 5) + dummy)*10)
  > 
  > pdata$x3[c(1,5,10)] <- NA
  > pdata$y[15] <- 1     # make things unbalanced
  > 
  > set.seed(10)
  > pfit <- rpart(y ~ x1 + x2 + x3, pdata,
  +               cp=0, xval=0, minsplit=5, maxdepth=1,
  +               parms=list(prior=c(.2, .3, .5), 
  +                          loss =matrix(c(0,2,2,2,0,6,1,1,0), 3,3,byrow=T)))
  > 
  > #
  > # See section 12.1 of the report for these numbers
  > #
  > 
  > ntot   <- c(6,5,4)
  > phat   <- c(6,5,4)/15   # observed class probabilities
  > prior  <- c(.2, .3, .5) # priors
  > aprior <- c(4,12,5)/21  # altered priors
  > lmat   <- matrix(c(0,1,2, 2,0,1, 2,6,0), ncol=3)  #loss matrix
  > 
  > gini <- function(p) 1-sum(p^2)
  > loss <- function(n, class)  sum(n * lmat[,class])
  > 
  > phat <- function(n, ntot=c(6,5,4), prior=c(.2, .3, .5)) {
  +     n*prior/ntot
  +     }
  > 
  > # Are the losses correct?
  > # Class counts for the two children are (4,4,0) and (2,1,4), when
  > #   using surrogates
  > aeq(pfit$frame$dev/15, c(loss(prior,2), loss(phat(c(4,4,0)),2), 
  +                                  loss(phat(c(2,1,4)),3)))
  [1] TRUE
  > # Node probabilities?
  > aeq(pfit$frame$yval2[,8] , 
  +     c(1, sum(phat(c(4,4,0))), sum(phat(c(2,1,4)))))
  [1] TRUE
  > 
  > aeq(pfit$frame$yval2[,5:7] , rbind(prior, 
  +                                    phat(c(4,4,0))/ sum(phat(c(4,4,0))), 
  +                                    phat(c(2,1,4))/ sum(phat(c(2,1,4)))))
  [1] TRUE
  > 
  > # Now the node and class probs, under altered priors
  > phat2 <- function(n, ntot=c(6,5,4), prior=aprior) {
  +     n*prior/ntot
  +     }
  > 
  > # Use these to create the gini losses, base data, and for the best
  > #  splits on variables 1, 2, 3
  > gfun <- function(n) {  #The gini loss for a node, given the counts
  +     temp <- phat2(n)
  +     sum(temp) * gini(temp/sum(temp))
  +     }
  >            
  > # These are in order x3, x2, x1 (best split to worst)
  > # Note that for x3, missing values cause the "parent" to be viewed as
  > #   having 12 obs instead of 15.
  > # Each line is gini(parent) - gini(children)
  > aeq(pfit$splits[1:3, 3],
  +     15* c(gfun(c(4,4,4)) - (gfun(c(3,4,0)) +  gfun(c(1,0,4))),
  +           gfun(c(6,5,4)) - (gfun(c(6,5,2)) +  gfun(c(0,0,2))),
  +           gfun(c(6,5,4)) - (gfun(c(4,4,4)) +  gfun(c(2,1,0)))))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.06    0.21 
Running the tests in 'tests/treble.R' failed.
Complete output:
  > #
  > # Simplest weight test: treble the weights
  > #
  > #  By using the unshrunken estimates the weights will nearly cancel
  > #   out:  frame$wt, frame$dev, frame$yval2, and improvement will all
  > #   be threefold larger, other things will be the same.
  > # The improvement is the splits matrix, column 3, rows with n>0.  Other
  > #   rows are surrogate splits.
  > library(rpart)
  > require(survival)
  Loading required package: survival
  > set.seed(10)
  > 
  > tempc <- rpart.control(maxsurrogate=0, cp=0, xval=0)
  > fit1 <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  +                 stagec, control=tempc,
  +                 method='poisson', parms=list(shrink=0))
  > wts <- rep(3, nrow(stagec))
  > fit1b <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  +                 stagec, control= tempc, parms=list(shrink=0),
  +                 method='poisson', weights=wts)
  > fit1b$frame$wt   <- fit1b$frame$wt/3
  > fit1b$frame$dev  <- fit1b$frame$dev/3
  > fit1b$frame$yval2[,2] <- fit1b$frame$yval2[,2]/3
  > fit1b$splits[,3] <- fit1b$splits[,3]/3
  > zz <- match(c("call", "variable.importance"), names(fit1))
  > all.equal(fit1[-zz], fit1b[-zz])   #all but the "call" and importance
  [1] TRUE
  > all.equal(fit1b$variable.importance/fit1$variable.importance, rep(3,4),
  +           check.attributes = FALSE)
  [1] TRUE
  > 
  > #
  > # Compare a pair of multiply weighted fits
  > #  In this one, the lengths of where and y won't match
  > # I have to set minsplit to the smallest possible, because otherwise
  > #  the replicated data set will sometimes have enough "n" to split, but
  > #  the weighted one won't.  Use of CP keeps the degenerate splits
  > #  (n=2, several covariates with exactly the same improvement) at bay.
  > # For larger trees, the weighted split will sometimes have fewer
  > #  surrogates, because of the "at least two obs" rule.
  > #
  > # Create a reproducable psuedo random order using the logisic attractor
  > pseudo <- double(nrow(stagec))
  > pseudo[1] <- pi/4
  > for (i in 2:nrow(stagec)) pseudo[i] <- 4*pseudo[i-1]*(1 - pseudo[i-1])
  > 
  > wts <- rep(1:5, length=nrow(stagec))
  > temp <- rep(1:nrow(stagec), wts)             #row replicates
  > xgrp <- rep(1:10, length=146)[order(pseudo)]
  > xgrp2<- rep(xgrp, wts)
  > #  Direct: replicate rows in the data set, and use unweighted
  > fit2 <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  + 	      control=rpart.control(minsplit=2, xval=xgrp2, cp=.025),
  + 	      data=stagec[temp,], method='poisson')
  > 
  > #  Weighted
  > fit2b<- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  + 	      control=rpart.control(minsplit=2, xval=xgrp, cp=.025),
  + 	      data=stagec, method='poisson', weight=wts)
  > 
  > all.equal(fit2$frame[-2],  fit2b$frame[-2])  # the "n" component won't match
  [1] TRUE
  > all.equal(fit2$cptable,    fit2b$cptable)
  [1] TRUE
  > #all.equal(fit2$splits[,-1],fit2b$splits[,-1]) #fails
  > toss <- c(49, 64)
  > all.equal(fit2$splits[-toss,-1],fit2b$splits[-toss,-1]) #ok
  [1] TRUE
  > all.equal(fit2$csplit,    fit2b$csplit)
  [1] TRUE
  > # Line 49 is a surrogate split in a group whose 2 smallest ages are
  > #  47 and 48.  The weighted fit won't split there because it wants to
  > #  send at least 2 obs to the left; the replicate fit thinks that there
  > #  are several 47's.
  > 
  > 
  > 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.17    1.53 
Running the tests in 'tests/treble2.R' failed.
Complete output:
  > #
  > # Test weights in a regression problem
  > #
  > library(rpart)
  > set.seed(10)
  > 
  > mystate <- data.frame(state.x77, region=factor(state.region))
  > names(mystate) <- c("population","income" , "illiteracy","life" ,
  +        "murder", "hs.grad", "frost",     "area",      "region")
  > 
  > xgrp <- rep(1:10,5)
  > fit4 <- rpart(income ~ population + region + illiteracy +life + murder +
  +                         hs.grad + frost , mystate,
  +                    control=rpart.control(minsplit=10, xval=xgrp))
  > wts <- rep(3, nrow(mystate))
  > fit4b <-  rpart(income ~ population + region + illiteracy +life + murder +
  +                         hs.grad + frost , mystate,
  +                    control=rpart.control(minsplit=10, xval=xgrp), weights=wts)
  > fit4b$frame$wt   <- fit4b$frame$wt/3
  > fit4b$frame$dev  <- fit4b$frame$dev/3
  > fit4b$cptable[,5] <- fit4b$cptable[,5] * sqrt(3)
  > temp <- c('frame', 'where', 'splits', 'csplit', 'cptable')
  > all.equal(fit4[temp], fit4b[temp])  
  [1] TRUE
  > 
  > 
  > # Next is a very simple case, but worth keeping
  > dummy <- data.frame(y=1:10, x1=c(10:4, 1:3), x2=c(1,3,5,7,9,2,4,6,8,0))
  > 
  > xx1 <- rpart(y ~ x1 + x2, dummy, minsplit=4, xval=0)
  > xx2 <- rpart(y ~ x1 + x2, dummy, weights=rep(2,10), minsplit=4, xval=0)
  > 
  > all.equal(xx1$frame$dev, c(82.5, 10, 2, .5, 10, .5, 2))
  [1] TRUE
  > all.equal(xx2$frame$dev, c(82.5, 10, 2, .5, 10, .5, 2)*2)
  [1] TRUE
  > 
  > # Now for a set of non-equal weights
  > #  We need to set maxcompete=3 because there just happens to be, in one
  > #  of the lower nodes, an exact tie between variables "life" and "murder".
  > #  Round off error causes fit5 to choose one and fit5b the other.
  > # Later -- cut it back to maxdepth=3 for the same reason (a tie).
  > #
  > nn <- nrow(mystate)
  > wts <- rep(1:5, length=nn)
  > temp <- rep(1:nn, wts)             #row replicates
  > xgrp <- rep(1:10, length=nn)
  > xgrp2<- rep(xgrp, wts)
  > tempc <- rpart.control(minsplit=2, xval=xgrp2, maxsurrogate=0,
  + 		       maxcompete=3, maxdepth=3)
  > #  Direct: replicate rows in the data set, and use unweighted
  > fit5 <-  rpart(income ~ population + region + illiteracy +life + murder +
  +                         hs.grad + frost , data=mystate[temp,], control=tempc)
  > #  Weighted
  > tempc <- rpart.control(minsplit=2, xval=xgrp, maxsurrogate=0,
  + 		       maxcompete=3, maxdepth=3)
  > fit5b <-  rpart(income ~ population + region + illiteracy +life + murder +
  +                         hs.grad + frost , data=mystate, control=tempc,
  +                         weights=wts)
  > all.equal(fit5$frame[-2],  fit5b$frame[-2])  # the "n" component won't match
  [1] TRUE
  > all.equal(fit5$cptable,    fit5b$cptable)
  [1] TRUE
  > all.equal(fit5$splits[,-1],fit5b$splits[,-1]) 
  [1] TRUE
  > all.equal(fit5$csplit,    fit5b$csplit)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.15    0.12    0.26 
Running the tests in 'tests/treble3.R' failed.
Complete output:
  > #
  > # The treble test for classification trees
  > #  
  > #
  > library(rpart)
  > set.seed(10)
  > 
  > xgrp <- rep(1:10,length=nrow(cu.summary))
  > carfit <- rpart(Country ~ Reliability + Price + Mileage + Type,
  + 		 method='class', data=cu.summary, 
  + 		 control=rpart.control(xval=xgrp))
  > 
  > carfit2 <- rpart(Country ~ Reliability + Price + Mileage + Type,
  + 		 method='class', data=cu.summary, 
  + 		 weight=rep(3,nrow(cu.summary)),
  + 		 control=rpart.control(xval=xgrp))
  > 
  > all.equal(carfit$frame$wt,    carfit2$frame$wt/3)
  [1] TRUE
  > all.equal(carfit$frame$dev,   carfit2$frame$dev/3)
  [1] TRUE
  > all.equal(carfit$frame[,5:7], carfit2$frame[,5:7])
  [1] TRUE
  > all.equal(carfit$frame$yval2[,12:21], carfit2$frame$yval2[,12:21])
  [1] TRUE
  > all.equal(carfit[c('where', 'csplit')],
  + 	  carfit2[c('where', 'csplit')])
  [1] TRUE
  > xx <- carfit2$splits
  > xx[,'improve'] <- xx[,'improve'] / ifelse(xx[,5]> 0,1,3) # surrogate?
  > all.equal(xx, carfit$splits)
  [1] TRUE
  > all.equal(as.vector(carfit$cptable), 
  + 	  as.vector(carfit2$cptable%*% diag(c(1,1,1,1,sqrt(3)))))
  [1] TRUE
  > 
  > summary(carfit2)
  Call:
  rpart(formula = Country ~ Reliability + Price + Mileage + Type, 
      data = cu.summary, weights = rep(3, nrow(cu.summary)), method = "class", 
      control = rpart.control(xval = xgrp))
    n= 117 
  
            CP nsplit rel error    xerror       xstd
  1 0.23529412      0 1.0000000 1.0000000 0.04530958
  2 0.02941176      1 0.7647059 0.8088235 0.04583685
  3 0.02205882      3 0.7058824 0.9117647 0.04583685
  4 0.01000000      5 0.6617647 0.9558824 0.04563477
  
  Variable importance
  Reliability       Price        Type 
           49          30          21 
  
  Node number 1: 117 observations,    complexity param=0.2352941
    predicted class=USA        expected loss=0.5811966  P(node) =1
      class counts:     3     3     6    33    93    27    15     9    15   147
     probabilities: 0.009 0.009 0.017 0.094 0.265 0.077 0.043 0.026 0.043 0.419 
    left son=2 (33 obs) right son=3 (84 obs)
    Primary splits:
        Reliability splits as  RRRLL,       improve=44.777630, (32 missing)
        Type        splits as  LRRLRR,      improve=11.246350, (0 missing)
        Price       < 6923    to the left,  improve=10.448720, (0 missing)
        Mileage     < 22.5    to the right, improve= 6.084489, (57 missing)
    Surrogate splits:
        Type  splits as  RRRLRR,      agree=0.706, adj=0.138, (32 split)
        Price < 6665    to the left,  agree=0.682, adj=0.069, (0 split)
  
  Node number 2: 33 observations,    complexity param=0.02205882
    predicted class=Japan      expected loss=0.4545455  P(node) =0.2820513
      class counts:     3     0     0     3    54    27     3     3     0     6
     probabilities: 0.030 0.000 0.000 0.030 0.545 0.273 0.030 0.030 0.000 0.061 
    left son=4 (9 obs) right son=5 (24 obs)
    Primary splits:
        Price       < 13202   to the right, improve=7.621212, (0 missing)
        Type        splits as  R-LRLL,      improve=5.621212, (0 missing)
        Reliability splits as  RRRRL,       improve=5.406404, (4 missing)
        Mileage     < 23.5    to the left,  improve=3.915966, (9 missing)
    Surrogate splits:
        Type splits as  R-LRLL, agree=0.879, adj=0.556, (0 split)
  
  Node number 3: 84 observations,    complexity param=0.02941176
    predicted class=USA        expected loss=0.4404762  P(node) =0.7179487
      class counts:     0     3     6    30    39     0    12     6    15   141
     probabilities: 0.000 0.012 0.024 0.119 0.155 0.000 0.048 0.024 0.060 0.560 
    left son=6 (30 obs) right son=7 (54 obs)
    Primary splits:
        Price       < 16668.5 to the right, improve=10.403170, (0 missing)
        Type        splits as  RRRLRR,      improve= 6.025974, (0 missing)
        Mileage     < 22.5    to the right, improve= 3.132611, (48 missing)
        Reliability splits as  RLLLL,       improve= 2.195489, (28 missing)
    Surrogate splits:
        Type splits as  RLLRRR, agree=0.714, adj=0.2, (0 split)
  
  Node number 4: 9 observations
    predicted class=Japan      expected loss=0.1111111  P(node) =0.07692308
      class counts:     0     0     0     3    24     0     0     0     0     0
     probabilities: 0.000 0.000 0.000 0.111 0.889 0.000 0.000 0.000 0.000 0.000 
  
  Node number 5: 24 observations,    complexity param=0.02205882
    predicted class=Japan      expected loss=0.5833333  P(node) =0.2051282
      class counts:     3     0     0     0    30    27     3     3     0     6
     probabilities: 0.042 0.000 0.000 0.000 0.417 0.375 0.042 0.042 0.000 0.083 
    left son=10 (7 obs) right son=11 (17 obs)
    Primary splits:
        Price   < 7038    to the left,  improve=4.836134, (0 missing)
        Type    splits as  R--LR-,      improve=1.966667, (0 missing)
        Mileage < 30.5    to the right, improve=1.716667, (6 missing)
  
  Node number 6: 30 observations,    complexity param=0.02941176
    predicted class=USA        expected loss=0.6333333  P(node) =0.2564103
      class counts:     0     3     3    24    12     0     0     0    15    33
     probabilities: 0.000 0.033 0.033 0.267 0.133 0.000 0.000 0.000 0.167 0.367 
    left son=12 (10 obs) right son=13 (20 obs)
    Primary splits:
        Type        splits as  LRR-L-,      improve=6.6000000, (0 missing)
        Price       < 30511   to the right, improve=5.3739130, (0 missing)
        Reliability splits as  RRLLL,       improve=0.3333333, (12 missing)
    Surrogate splits:
        Price < 19422.5 to the left,  agree=0.767, adj=0.3, (0 split)
  
  Node number 7: 54 observations
    predicted class=USA        expected loss=0.3333333  P(node) =0.4615385
      class counts:     0     0     3     6    27     0    12     6     0   108
     probabilities: 0.000 0.000 0.019 0.037 0.167 0.000 0.074 0.037 0.000 0.667 
  
  Node number 10: 7 observations
    predicted class=Japan      expected loss=0.2857143  P(node) =0.05982906
      class counts:     0     0     0     0    15     3     3     0     0     0
     probabilities: 0.000 0.000 0.000 0.000 0.714 0.143 0.143 0.000 0.000 0.000 
  
  Node number 11: 17 observations
    predicted class=Japan/USA  expected loss=0.5294118  P(node) =0.1452991
      class counts:     3     0     0     0    15    24     0     3     0     6
     probabilities: 0.059 0.000 0.000 0.000 0.294 0.471 0.000 0.059 0.000 0.118 
  
  Node number 12: 10 observations
    predicted class=Germany    expected loss=0.5  P(node) =0.08547009
      class counts:     0     0     0    15     3     0     0     0     9     3
     probabilities: 0.000 0.000 0.000 0.500 0.100 0.000 0.000 0.000 0.300 0.100 
  
  Node number 13: 20 observations
    predicted class=USA        expected loss=0.5  P(node) =0.1709402
      class counts:     0     3     3     9     9     0     0     0     6    30
     probabilities: 0.000 0.050 0.050 0.150 0.150 0.000 0.000 0.000 0.100 0.500 
  
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.09    0.25 
Running the tests in 'tests/usersplits.R' failed.
Complete output:
  > # Any necessary setup
  > library(rpart)
  > options(na.action="na.omit")
  > options(digits=4) # to match earlier output
  > set.seed(1234)
  > 
  > mystate <- data.frame(state.x77, region=factor(state.region))
  > names(mystate) <- c("population","income" , "illiteracy","life" ,
  +        "murder", "hs.grad", "frost",     "area",      "region")
  > #
  > # Test out the "user mode" functions, with an anova variant
  > #
  > 
  > # The 'evaluation' function.  Called once per node.
  > #  Produce a label (1 or more elements long) for labeling each node,
  > #  and a deviance.  The latter is
  > #	- of length 1
  > #       - equal to 0 if the node is "pure" in some sense (unsplittable)
  > #       - does not need to be a deviance: any measure that gets larger
  > #            as the node is less acceptable is fine.
  > #       - the measure underlies cost-complexity pruning, however
  > temp1 <- function(y, wt, parms) {
  +     wmean <- sum(y*wt)/sum(wt)
  +     rss <- sum(wt*(y-wmean)^2)
  +     list(label= wmean, deviance=rss)
  +     }
  > 
  > # The split function, where most of the work occurs.
  > #   Called once per split variable per node.
  > # If continuous=T
  > #   The actual x variable is ordered
  > #   y is supplied in the sort order of x, with no missings,
  > #   return two vectors of length (n-1):
  > #      goodness = goodness of the split, larger numbers are better.
  > #                 0 = couldn't find any worthwhile split
  > #        the ith value of goodness evaluates splitting obs 1:i vs (i+1):n
  > #      direction= -1 = send "y< cutpoint" to the left side of the tree
  > #                  1 = send "y< cutpoint" to the right
  > #         this is not a big deal, but making larger "mean y's" move towards
  > #         the right of the tree, as we do here, seems to make it easier to
  > #         read
  > # If continuos=F, x is a set of integers defining the groups for an
  > #   unordered predictor.  In this case:
  > #       direction = a vector of length m= "# groups".  It asserts that the
  > #           best split can be found by lining the groups up in this order
  > #           and going from left to right, so that only m-1 splits need to
  > #           be evaluated rather than 2^(m-1)
  > #       goodness = m-1 values, as before.
  > #
  > # The reason for returning a vector of goodness is that the C routine
  > #   enforces the "minbucket" constraint. It selects the best return value
  > #   that is not too close to an edge.
  > temp2 <- function(y, wt, x, parms, continuous) {
  +     # Center y
  +     n <- length(y)
  +     y <- y- sum(y*wt)/sum(wt)
  + 
  +     if (continuous) {
  + 	# continuous x variable
  + 	temp <- cumsum(y*wt)[-n]
  + 
  + 	left.wt  <- cumsum(wt)[-n]
  + 	right.wt <- sum(wt) - left.wt
  + 	lmean <- temp/left.wt
  + 	rmean <- -temp/right.wt
  + 	goodness <- (left.wt*lmean^2 + right.wt*rmean^2)/sum(wt*y^2)
  + 	list(goodness= goodness, direction=sign(lmean))
  + 	}
  +     else {
  + 	# Categorical X variable
  + 	ux <- sort(unique(x))
  + 	wtsum <- tapply(wt, x, sum)
  + 	ysum  <- tapply(y*wt, x, sum)
  + 	means <- ysum/wtsum
  + 
  + 	# For anova splits, we can order the categories by their means
  + 	#  then use the same code as for a non-categorical
  + 	ord <- order(means)
  + 	n <- length(ord)
  + 	temp <- cumsum(ysum[ord])[-n]
  + 	left.wt  <- cumsum(wtsum[ord])[-n]
  + 	right.wt <- sum(wt) - left.wt
  + 	lmean <- temp/left.wt
  + 	rmean <- -temp/right.wt
  + 	list(goodness= (left.wt*lmean^2 + right.wt*rmean^2)/sum(wt*y^2),
  + 	     direction = ux[ord])
  + 	}
  +     }
  > 
  > # The init function:
  > #   fix up y to deal with offsets
  > #   return a dummy parms list
  > #   numresp is the number of values produced by the eval routine's "label"
  > #   numy is the number of columns for y
  > #   summary is a function used to print one line in summary.rpart
  > # In general, this function would also check for bad data, see rpart.poisson
  > #   for instace.
  > temp3 <- function(y, offset, parms, wt) {
  +     if (!is.null(offset)) y <- y-offset
  +     list(y=y, parms=0, numresp=1, numy=1,
  + 	      summary= function(yval, dev, wt, ylevel, digits ) {
  + 		  paste("  mean=", format(signif(yval, digits)),
  + 			", MSE=" , format(signif(dev/wt, digits)),
  + 			sep='')
  + 	     })
  +     }
  > 
  > 
  > alist <- list(eval=temp1, split=temp2, init=temp3)
  > 
  > fit1 <- rpart(income ~population +illiteracy  + murder + hs.grad + region,
  + 	     mystate, control=rpart.control(minsplit=10, xval=0),
  + 	     method=alist)
  > 
  > fit2 <- rpart(income ~population +illiteracy + murder + hs.grad + region,
  + 	     mystate, control=rpart.control(minsplit=10, xval=0),
  + 	      method='anova')
  > 
  > # Other than their call statement, and a longer "functions" component in
  > #  fit1, fit1 and fit2 should be identical.
  > all.equal(fit1$frame, fit2$frame)
  [1] TRUE
  > all.equal(fit1$splits, fit2$splits)
  [1] TRUE
  > all.equal(fit1$csplit, fit2$csplit)
  [1] TRUE
  > all.equal(fit1$where, fit2$where)
  [1] TRUE
  > all.equal(fit1$cptable, fit2$cptable)
  [1] TRUE
  > 
  > # Now try xpred on it
  > xvtemp <- rep(1:5, length=50)
  > xp1 <- xpred.rpart(fit1, xval=xvtemp)
  > xp2 <- xpred.rpart(fit2, xval=xvtemp)
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > aeq(xp1, xp2)
  [1] TRUE
  > 
  > fit3 <- rpart(income ~population +illiteracy + murder + hs.grad + region,
  + 	     mystate, control=rpart.control(minsplit=10, xval=xvtemp),
  + 	      method='anova')
  > zz <- apply((mystate$income - xp1)^2,2, sum)
  > aeq(zz/fit1$frame$dev[1], fit3$cptable[,4])  #reproduce xerror
  [1] TRUE
  > 
  > zz2 <- sweep((mystate$income-xp1)^2,2, zz/nrow(xp1))
  > zz2 <- sqrt(apply(zz2^2, 2, sum))/ fit1$frame$dev[1]
  > aeq(zz2, fit3$cptable[,5])          #reproduce se(xerror)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.34    0.04    0.34 
Running the tests in 'tests/xpred2.R' failed.
Complete output:
  > #
  > # Test out the "return.all" argument of xpred
  > #  The data set has the virtue of continuous, categorical, and missings
  > #
  > library(rpart)
  > require(survival)
  Loading required package: survival
  > set.seed(10)
  > 
  > fit1 <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  +                 stagec,  method='poisson')
  > 
  > xgrp <- rep(1:3, length=nrow(stagec))  # explicitly set the xval groups
  > 
  > xfit1 <- xpred.rpart(fit1, xval=xgrp, return.all=T)
  > xfit2 <- array(0, dim=dim(xfit1))
  > cplist <- as.numeric(dimnames(xfit1)[[2]])
  > 
  > for (i in 1:3) {
  +     tfit <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2+grade+gleason +ploidy,
  +                 stagec,  method='poisson', subset=(xgrp !=i))
  +     # xvals are actually done on the absolute risk (node's risk /n), not on
  +     #   the rescaled risk ((node risk)/ (top node risk)) which is the basis
  +     #   for the printed CP.  To get the right answer we need to rescale.
  +     cp2 <- cplist * (fit1$frame$dev[1] / fit1$frame$n[1]) / 
  +                     (tfit$frame$dev[1] / tfit$frame$n[1])
  + 
  +     for (j in 1:length(cp2)) {
  +         tfit2 <- prune(tfit, cp=cp2[j])
  +         temp <- predict(tfit2, newdata=stagec[xgrp==i,], type='matrix')
  +         xfit2[xgrp==i, j,] <- temp
  +         }
  +     }
  > 
  > all.equal(xfit1, xfit2, check.attributes=FALSE)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.40    0.20    1.59