CI after logistic regression in bootstrap

General
1

Hi, the following codes provided me results of bootstrapped logistic regression, but boot.ci() only provide one CI, which I have no idea what it it. I wonder if there is a way to obtain the CI for each factor in the model. Any help is appreciated.
logistic_coef <-function(d, i){
d <- d[i,]
fit <- glm(bpd ~ birthwt + gestage + toxemia,
data = d,
family = "binomial")
return(coef(fit))
}

#apply the boot() function
#sytem.time estimate time used
set.seed(77)
system.time(boot_logistic <- boot(data = bpd,
statistic = logistic_coef,
R=5000))
#view results
boot_logistic
boot.ci(boot_logistic, type = "perc")
Call:
boot(data = bpd, statistic = logistic_coef, R = 5000)

Bootstrap Statistics :
original bias std. error
t1* 13.936082577 0.4182376970 2.7596107013
t2* -0.002643578 -0.0001080289 0.0009459613
t3* -0.388535681 -0.0107899365 0.1042887437
t4* -1.343786469 -0.1162307066 1.0181248939

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates

CALL :
boot.ci(boot.out = boot_logistic, type = "perc")

Intervals :
Level Percentile
95% ( 9.32, 20.01 )
Calculations and Intervals on Original Scale

2

Hi @Jason1,

You are encouraged to make your question reproducible (see this article: FAQ: How to do a minimal reproducible example ( reprex ) for beginners). If you do, it will be easier for the community here to help you out. In your case, I believe you only need to share a sample dataset so we can run your code.

1 Like
3

Hi Gueyenono,
I've pasted the first 50 observations of the data below. your help is much appreciated.
|bpd|birthwt|gestage|toxemia|

1 850 27 0
0 1500 33 0
1 1360 32 0
0 960 35 1
0 1560 33 0
0 1120 29 0
1 810 28 0
0 1620 32 0
1 1000 30 0
1 700 26 0
0 1330 31 0
0 1410 31 0
1 1520 31 0
1 910 29 0
0 1650 33 0
0 1460 32 0
1 1000 28 0
1 710 25 0
1 1220 28 0
1 820 28 0
0 1060 29 0
0 1240 30 0
0 1440 33 0
0 1470 32 0
0 790 29 0
1 910 29 0
0 1020 29 0
0 1340 28 0
0 1140 34 0
0 1450 34 0
1 1670 30 0
0 1690 33 0
0 950 31 0
1 810 29 0
1 730 27 1
1 730 26 0
1 650 30 0
0 900 26 0
0 1690 32 0
0 1200 28 0
1 700 26 0
1 750 29 0
1 1000 31 0
0 1500 31 0
1 1310 31 0
1 1310 31 0
1 990 29 0
1 910 28 0
0 1500 33 0
1 1030 27 0
4

@Jason1 it's not easy to use your data when it is pasted in this format. If you imported your data in mydata for example, run dput(my data) and post your result here.

5

Hi Gueyenono,

here are the data as suggested. thanks.
dput(bpdsample)
structure(list(bpd = structure(c(1L, 1L, 2L, 2L, 1L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L,
1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L), .Label = c("0",
"1"), class = "factor"), birthwt = c(870, 1340, 850, 1420, 1250,
930, 1310, 1330, 1470, 1150, 1570, 1560, 700, 820, 1680, 1180,
520, 1650, 1170, 650, 1260, 1180, 1040, 1690, 1380, 1130, 730,
1400, 1670, 1530, 1550, 1350, 1520, 820, 1200, 820, 820, 1500,
1680, 990, 700, 1050, 1000, 870, 1320, 1380, 1550, 890, 1410,
1050), gestage = c(29, 30, 30, 30, 32, 28, 31, 30, 32, 32, 33,
36, 26, 28, 30, 32, 28, 33, 29, 26, 31, 30, 29, 33, 32, 30, 27,
29, 36, 31, 33, 31, 33, 29, 30, 26, 28, 30, 32, 30, 26, 28, 31,
28, 31, 33, 32, 27, 31, 29), toxemia = structure(c(1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L), .Label = c("0",
"1"), class = "factor")), row.names = c("97", "155", "166", "99",
"148", "57", "45", "154", "24", "184", "76", "211", "10", "181",
"205", "199", "165", "158", "74", "81", "91", "70", "54", "32",
"147", "69", "35", "138", "84", "130", "159", "153", "168", "109",
"80", "163", "20", "95", "216", "221", "41", "144", "43", "161",
"92", "220", "68", "63", "12", "136"), class = "data.frame")

6

There is no problem with your results. They already contain bootstrap statistics for each of the coefficients in your logistic model.: t1* is the intercept, t2* is for birthwt, t3* is for gestage and t4* is for toxemial.

7

I got that, but my question is how to get the confidence interval for each of them. when using bent(), it provides only 1 CI, which I guess it is for the intercept. Do you know how to get the rest of them? Thanks.

8

Sorry, I mean bent.ci().

9

Hi, I think I found the solution, this seems to be what you need to get the CI of each variable
boot.ci(boot_logistic, type = "perc",index = 1)

10

You mentioned bent.ci() earlier and I was not sure what function you were referring to. It turns out that you use boot.ci() once again in your solution. What is the difference between your solution now and what you got from the very beginning?

11

you can replace the # in "index =#" to get each CI.

12

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.