【问题标题】:bootstrap to get confidence intervals using R使用 R 引导获取置信区间
【发布时间】:2017-07-15 18:10:23
【问题描述】:

我正在尝试使用非参数引导程序引导可靠性估计 我已经编写了下面的代码,其中创建了一个模型,然后引导 1000 次以获得两个可靠性统计 Alpha 和 Omega 我能够获得具有置信区间的第一个构造的 Alpha 和 Omega:visual =~ x1 + x2 + x3 但看不到其他构造 textualspeed 的访问方法 当我运行引导功能时,我会看到所有结果

# 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

以下是我公认的粗制滥造的尝试。谁能帮忙

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)]

感谢您的帮助

【问题讨论】:

  • 基本调试...在新会话中执行此操作...并阅读错误消息。我遇到的第一个错误是: ERsum(beta[i, ], tau.found.sym.optim, m + 1, m + n) 中的错误:dims [product 666] do not match the length of object [999] In另外:警告消息:在 y - X %*% beta 中:较长的对象长度不是较短对象长度的倍数
  • 对不起,是的,我现在看到了。我会马上更新代码

标签: r statistics-bootstrap


【解决方案1】:

您首先需要查看从 reliability 返回的完整数据集:

> 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

那么你需要看看你的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"

....所以results$t0 包含所有三个模型参数估计值。

【讨论】:

  • 嗨@42-,谢谢你的回复。我能够提取omegaalpha 点统计信息,但我也想获得当我为所有三个构造运行代码boot.ci(results, type="bca", index = 1) 时我无法看到的置信区间; visualtextualspeed
  • 阅读?boot.ci 的帮助页面并特别注意“index”参数。使用默认值,您只能获得第一个参数的 CI 估计值。
猜你喜欢
  • 2016-11-28
  • 1970-01-01
  • 2019-05-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2018-07-29
  • 1970-01-01
相关资源
最近更新 更多