[R] boot statictic fn for dual estimation of 2 stats?  
    Olivia Lau 
    olau at fas.harvard.edu
       
    Sat Oct 11 21:49:10 CEST 2003
    
    
  
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)
    
    
More information about the R-help
mailing list