bootstrap to get confidence intervals using R

Question

Im trying to bootstrap reliability estimates using non parametric bootstrap I have written the code below where a model is created and then bootstrapped 1000 times in order to get two reliability statistics Alpha and Omega I am able to get the Alpha and the Omega for the first construct with confidence intervals: visual =~ x1 + x2 + x3 but see no way of accessing it for the other constructs textual and speed When i run the boot function, i see results for all of them

# bootstrapping with 1000 replications
results <- boot(data=data, statistic=reliability, R=500, formula=HS.model,parallel = 'snow')

> results$t0

        visual   textual     speed     total
alpha  0.6261171 0.8827069 0.6884550 0.7604886
omega  0.6253180 0.8851754 0.6877600 0.8453351
omega2 0.6253180 0.8851754 0.6877600 0.8453351
omega3 0.6120052 0.8850608 0.6858417 0.8596204
avevar 0.3705589 0.7210163 0.4244883 0.5145874

Below is my admittedly shoddy attempt. Can anyone help

library(lavaan)
library(semTools)
library(boot)

data <- HolzingerSwineford1939

HS.model <- 'visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9 '

# function to reliability stats
reliability <- function(formula, data, indices) {
  data = data
  d <- data[indices,] # allows boot to select sample
  fit <- cfa(HS.model, data=d)
  semTools::reliability(fit)
}

# bootstrapping with 500 replications
results <- boot(data=data, statistic=reliability, R=500, formula=HS.model,parallel = 'snow')

# Get the confidence intervals
conf_interval_alpha <- boot.ci(results, type="bca", index = 1)

# Retrieve the Alpha and confidence intervals
alpha <- conf_interval_alpha$t0
alpha.ci <- conf_interval_alpha$bca[,c(4,5)]

# Retrieve the Omega and confidence intervals  
conf_interval_omega <- boot.ci(results, type="bca", index = 2)
omega <- conf_interval_omega$t0
omega.ci <- conf_interval_omega$bca[,c(4,5)]

Thank you for your help

Basic debugging ... do this in fresh session ...and READ error messages. The first error I encountered was: Error in ERsum(beta[i, ], tau.found.sym.optim, m + 1, m + n) : dims [product 666] do not match the length of object [999] In addition: Warning message: In y - X %*% beta : longer object length is not a multiple of shorter object length — IRTFM, Jul 15 '17 at 18:24
Sorry, yes, i see it now. I will update the code straight away — John Smith, Jul 15 '17 at 18:28

IRTFM · Answer 1 · 2017-07-16T18:12:52.883

You first need to look at what is returned from reliability with the full dataset:

> reliability(data=data)
          visual   textual     speed     total
alpha  0.6261171 0.8827069 0.6884550 0.7604886
omega  0.6253180 0.8851754 0.6877600 0.8453351
omega2 0.6253180 0.8851754 0.6877600 0.8453351
omega3 0.6120052 0.8850608 0.6858417 0.8596204
avevar 0.3705589 0.7210163 0.4244883 0.5145874

Then you need to see what is returned from your boot-call:

> str(results)
List of 11
 $ t0       : num [1:5, 1:4] 0.626 0.625 0.625 0.612 0.371 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "alpha" "omega" "omega2" "omega3" ...
  .. ..$ : chr [1:4] "visual" "textual" "speed" "total"
 $ t        : num [1:500, 1:20] 0.594 0.607 0.613 0.669 0.621 ...
 $ R        : num 500
 $ data     :'data.frame':  301 obs. of  15 variables:
  ..$ id    : int [1:301] 1 2 3 4 5 6 7 8 9 11 ...
  ..$ sex   : int [1:301] 1 2 2 1 2 2 1 2 2 2 ...
  ..$ ageyr : int [1:301] 13 13 13 13 12 14 12 12 13 12 ...
  ..$ agemo : int [1:301] 1 7 1 2 2 1 1 2 0 5 ...
  ..$ school: Factor w/ 2 levels "Grant-White",..: 2 2 2 2 2 2 2 2 2 2 ...
  ..$ grade : int [1:301] 7 7 7 7 7 7 7 7 7 7 ...
  ..$ x1    : num [1:301] 3.33 5.33 4.5 5.33 4.83 ...
  ..$ x2    : num [1:301] 7.75 5.25 5.25 7.75 4.75 5 6 6.25 5.75 5.25 ...
  ..$ x3    : num [1:301] 0.375 2.125 1.875 3 0.875 ...
  ..$ x4    : num [1:301] 2.33 1.67 1 2.67 2.67 ...
  ..$ x5    : num [1:301] 5.75 3 1.75 4.5 4 3 6 4.25 5.75 5 ...
  ..$ x6    : num [1:301] 1.286 1.286 0.429 2.429 2.571 ...
  ..$ x7    : num [1:301] 3.39 3.78 3.26 3 3.7 ...
  ..$ x8    : num [1:301] 5.75 6.25 3.9 5.3 6.3 6.65 6.2 5.15 4.65 4.55 ...
  ..$ x9    : num [1:301] 6.36 7.92 4.42 4.86 5.92 ...
 $ seed     : int [1:626] 403 330 1346657232 1136157038 -874329217 857221657 1850455833 952027245 2020402269 -1198488986 ...
 $ statistic:function (formula, data, indices)  
  ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 16 6 1 16 1 1 6
  .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x7fa57918d430> 
 $ sim      : chr "ordinary"
 $ call     : language boot(data = data, statistic = reliability, R = 500, formula = HS.model, parallel = "snow")
 $ stype    : chr "i"
 $ strata   : num [1:301] 1 1 1 1 1 1 1 1 1 1 ...
 $ weights  : num [1:301] 0.00332 0.00332 0.00332 0.00332 0.00332 ...
 - attr(*, "class")= chr "boot"

.... so results$t0 has all three model parameter estimates in it.

Hi @42-, thank you for reply. I am able to pull the `omega` and the `alpha` point statistics but i also want to get the confidence intervals which im not able to see when i run the code `boot.ci(results, type="bca", index = 1)` for all three constructs; `visual`, `textual` and`speed` — John Smith, Jul 16 '17 at 15:02
Read the help page for `?boot.ci` and pay particular attention to the "index" parameter. With the default values you only get the CI estimate for the first parameter. — IRTFM, Jul 16 '17 at 18:10

bootstrap to get confidence intervals using R

1 Answers1