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
>