2
$\begingroup$

I am running a glm model using bootstrap, I can extract the coefficient mean and the confidence intervals for all the factors in my model.

But how can I get the pvalue from there?

Model:

rate.lm <- glm(one ~ two + factor(three), data = data)

lmfit <- function(data, indices) {
  fit <- glm(one ~ two + factor(three), data = data[indices, ])
  return(coef(fit))
}

results <- boot(data= data, statistic = lmfit, R = 1000, strata=data$two)

for (i in 2:length(rate.lm$coefficients)) {
      bci <- boot.ci(results, type = "basic", index = i)
      print(sprintf("%s,%.4f,%.4f,%.4f", names(rate.lm$coefficients)[i], results$t0[i],
                    bci$basic[4], bci$basic[5]))
}

Which gives me back the coeff mean and the 95% confidence interval. In that case I know the variable two is significant but I need a pvalue absolutely for the rest of my analysis and I am unsure on how to get it. Thanks for your help!

[1] "two,0.0101,-0.0022,0.0229"
[1] "factor(three)2,-0.2029,-0.2185,-0.1888"
[1] "factor(three)3,-0.1752,-0.1914,-0.1595"
[1] "factor(three)4,-0.2093,-0.2259,-0.1947"
[1] "factor(three)5,-0.1745,-0.1910,-0.1588"
[1] "factor(three)6,-0.1988,-0.2139,-0.1849"
$\endgroup$
1
  • $\begingroup$ Try running predFun1<-function(.)(coef(summary(.))[,"Pr(>|z|)"]) , and then replace coef(fit) with predFun1(fit) in your function $\endgroup$
    – Tom
    Oct 20, 2020 at 22:39

1 Answer 1

3
$\begingroup$

John Fox's Bootstrapping Regression Models is a great resource to get started as it provides lots of examples.

To answer your question directly, you'll simply need to access the t0 and t values provided by the boot object. To wit, boot$t0 is a numeric (of length $M+1$ where $M$ is the number of covariates) with the coefficients returned by running the regression with the original data and boot$t is a a matrix of dimension $R \times (M+1)$ holding the coefficients returned from a regression run on a bootstrapped sample of the data in each row.

Therefore, in your example, the one-sided p-values can be calculated pretty simply as follows:

R <- 1e3 # number of bootstraps

extrema <- apply(X=results$t, MARGIN=1, FUN=">", results$t0)

pvals <- rowSums(extrema)/(R+1)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.