[R] boot statictic fn for dual estimation of 2 stats?

[Roger D. Peng]

Occasionaly, the polr returns a zeta element of length equal to 2 rather 
than 3 (which I guess is what you're expecting).  You're bootfn function 
should always check to see that it's returning an object of the same 
length every time.

-roger

Olivia Lau wrote:
> Hi, 
> 
> I am trying to use boot() to refit an ordinal logit (polr in MASS) model.  
> (A very basic bootstrap which samples from the data frame without 
> replacement and updates the model.)
> 
> I need to extract two statistics per run (the coefficients and zeta) and I 
> tried concatenating them into a single vector after fitting, but I get the 
> following error: 
> 
> Error in "[<-"(*tmp*, r, , value = statistic(data, i[r, ], ...)) : 
> 	number of items to replace is not a multiple of replacement length
> 
> This error goes away if I just return the coefficients, but I need zeta in 
> order to calculate predicted probabilities for each outcome category (1:4 
> in this case).   Alternatively, if boot() could return a vector of 
> predicted probabilities for each categorical outcome, that would 
> work as well, but I don't have a clue as to how to start programming this. 
> 
> Thanks, Olivia Lau
> 
> ------------
> 
> Here's some sample data and code:  
> 
> cost <- c(4, 3, 2, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 4, 3, 1, 2, 2, 1, 3, 2, 
> 1, 4, 1, 1, 2, 2, 2, 1, 3, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 
> 2, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 3, 1, 3, 3, 1, 2, 1, 3, 2, 2, 2, 
> 3, 1, 2, 1, 2, 2, 1, 2)
> 
> mil <- c(1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
> 1, 0, 0, 0, 0, 1)
> 
> coop <-  c(4, 2, 1, 1, 1, 1, 2, 1, 2, 2, 1, 1, 3, 3, 3, 1, 4, 3, 1, 3, 4, 
> 1, 1, 1, 1, 3, 1, 1, 1, 4, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1, 
> 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 1, 3, 2, 3, 2, 3, 1, 1, 3, 2, 
> 2, 3, 2, 1, 4, 1, 3)
> 
> sanction <- data.frame(cost=cost, mil=mil, coop=coop)
> 
> est <- polr(as.factor(cost) ~ mil + coop, data = sanction)
> 
> bootfn <- function(data, i, object) {  
>   d <- data[i,]
>   fit <- update(object, data = d)
>   c(fit$coef, fit$zeta)
> }
> 
> res <- boot(sanction, bootfn, R = 10, object = est)
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>