[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