标签: statistics-bootstrap

我的目的是编写几个函数，旨在找到两个协方差矩阵之间的整体相似性，方法是将它们与随机向量相乘并关联响应向量，或者通过自举矩阵之一来获得可用于比较的相关系数分布。但在这两种情况下，我都得到了错误的结果。观察到的矩阵间相关性高达 0.93，但分布最多仅达到 0.2。这是函数的代码：

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE) …

如果我想使用boot()R's boot软件包中的函数来计算两个向量之间的Pearson相关系数的重要性,我应该这样做:

boot(re1, cor, R = 1000)

re1这两个观测向量的两列矩阵在哪里？我似乎无法正确得到这个因为cor这些向量0.8,但上面的函数返回-0.2为t0.

我今天在R代码中意识到了一种奇怪的行为.我尝试了一个包{boot.StepAIC},其中包含一个用于AIC逐步回归结果的bootstrap函数.但是我认为统计背景不是问题(我希望如此).
我可以使用R顶层的函数.这是我的示例代码.

require(MASS)
require(boot.StepAIC)

n<-100
x<-rnorm(n); y<-rnorm(n,sd=2); z<-rnorm(n,sd=3); res<-x+y+z+rnorm(n,sd=0.1)
dat.test<-as.data.frame(cbind(x,y,z,res))
form.1<-as.formula(res~x+y+z)
boot.stepAIC(lm(form.1, dat.test),dat.test) # should be OK - works at me

但是,我想把它包装在一个自己的函数中.我将数据和公式传递给该函数.但我在boot.stepAIC()中遇到错误:

100个引导样本中的模型拟合失败strsplit中的错误(nam.vars,":"):非字符参数

# custom function
fun.boot.lm.stepAIC<-function(dat,form) {
  if(!inherits(form, "formula")) stop("No formula given")
  fit.lm<-lm(formula=form,data=dat)
  return(boot.stepAIC(object=fit.lm,data=dat))
 }
fun.boot.lm.stepAIC(dat=dat.test,form=form.1)
# results in an error 

那错误在哪里？我想它必须与当地和全球环境有关,不是吗？

我非常感谢你有时间阅读这篇文章.

我有一个超大的30GB文件,包含600万条记录和3000条(主要是分类数据)列,采用csv格式.我想引导用于多项式回归的子样本,但即使我的机器中使用64GB RAM并且交换文件的两倍,它也很难实现,这个过程变得非常缓慢并且停止.

我正在考虑在R中生成子样本指标并使用sed或awk将它们输入系统命令,但不知道如何执行此操作.如果有人知道使用R命令干净的方法,我会非常感激.

一个问题是我需要选择子样本的完整观察,即我需要具有特定多项观察的所有行 - 它们从观察到观察的长度不同.我计划使用glmnet然后进行一些花哨的变换来得到多项式的近似值.另一点是我不知道如何选择样本大小以适应内存限制.

非常感谢你的想法.

R.version
platform       x86_64-pc-linux-gnu          
arch           x86_64                       
os             linux-gnu                    
system         x86_64, linux-gnu            
status                                      
major          2                            
minor          15.1                         
year           2012                         
month          06                           
day            22                           
svn rev        59600                        
language       R                            
version.string R version 2.15.1 (2012-06-22)
nickname       Roasted Marshmallows   

尤达

我想将bootstrap统计值(原始,偏差和错误)的值放到一个单独的列表中 - 但我无法弄清楚如何做到这一点.

这是一个例子:

> library(boot)
> set.seed(123)
> mean.fun <- function(data, idx) { mean(data[idx]) }
> data <- boot(data=rnorm(100), statistic=mean.fun, R=999)
> names(data)
 [1] "t0"        "t"         "R"         "data"     
 [5] "seed"      "statistic" "sim"       "call"     
 [9] "stype"     "strata"    "weights"  
> data
ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = rnorm(100), statistic = mean.fun, R = 999)


Bootstrap Statistics :
    original   bias    std. error
t1* 0.09040591 0.004751773  0.08823615

现在,而不是文本,我想要实际值.显然data$t0是"原始",但我不知道如何获得偏见和错误的值.

此外,由于输入函数名给你的代码,我输入boost的R和从源代码复制片段,并试图寻找它在我的本地研发设施.但找不到任何东西.为什么,R不应该从本地存储中获取源代码？

library(boot)
set.seed(1)
x=sample(0:1000,1000)
y=function(u,i) sum(x[i])
o=boot(x,y,1000)
theta1=NULL
theta1=cbind(theta1,o$t)
b=theta1[order(theta1)]
bp1=c(b[25], b[975])
ci=boot.ci(o,type="perc")

我使用两种方法来构建引导百分位数置信区间，但我得到了两个不同的答案。

bp1=c(b[25], b[975]) get (480474,517834)

同时ci=boot.ci(o,type="perc")得到 (480476, 517837 )

boot.ci 如何构建百分位区间？

我想得到自举的t值和lm的自举p值.我有以下代码(基本上是从一篇论文中复制而来).

# First of all you need the following packages
install.packages("car") 
install.packages("MASS")
install.packages("boot")
library("car")
library("MASS")
library("boot")

boot.function <- function(data, indices){
data <- data[indices,]
mod <- lm(prestige ~ income + education, data=data) # the liear model

# the first element of the following vector contains the t-value
# and the second element is the p-value
c(summary(mod)[["coefficients"]][2,3], summary(mod)[["coefficients"]][2,4])     
}

现在,我计算了bootstrapping模型,它给出了以下内容:

duncan.boot <- boot(Duncan, boot.function, 1999)
duncan.boot

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Duncan, statistic = boot.function, R = 1999)


Bootstrap Statistics : …

lme4我正在R 中使用混合模型：

full_mod3=lmer(logcptplus1 ~ logdepth*logcobb + (1|fyear) + (1 |flocation),
data=cpt, REML=TRUE)

概括：

Formula: logcptplus1 ~ logdepth * logcobb + (1 | fyear) + (1 | flocation)
   Data: cpt

REML criterion at convergence: 577.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7797 -0.5431  0.0248  0.6562  2.1733 

 Random effects:
 Groups    Name        Variance Std.Dev.
 fyear     (Intercept) 0.2254   0.4748  
 flocation (Intercept) 0.1557   0.3946  
 Residual              0.9663   0.9830  
Number of obs: 193, groups:  fyear, 16; flocation, 16

Fixed effects:
                  Estimate Std. Error t value …

我已经尝试过pip install scipy，一切看起来都很好，通过我打开文件的路径，找不到任何提及引导库的内容，尽管它位于他们的网站上：https ://docs.scipy.org/doc/scipy/reference/生成/scipy.stats.bootstrap.html

查看 Github https://github.com/scipy/scipy/blob/master/scipy/stats/_bootstrap.py后，我可以看到 5 天前有更新，尽管我三天前运行了代码，没有任何问题

我想使用新lme4包的新bootMer()功能(当前的开发人员版本).我是R的新手,不知道我应该为其FUN参数编写哪个函数.它说它需要一个数字向量,但我不知道该函数将执行什么.所以我有一个混合模型公式,它被转换为bootMer(),并且有许多重复.所以我不知道那个外部函数是做什么的？它应该是引导方法的模板吗？他的bootMer中是否已经实现了自举方法？那么为什么他们需要外部的"利益统计"呢？我应该使用哪种统计数据？

以下语法是否适用？R继续产生错误,产生FUN必须是数字向量.我不知道如何将估计与"适合"分开,甚至我应该首先做到这一点？我可以说我迷失了那种"有趣"的说法.另外我不知道我应该使用变量"Mixed5"传递混合模型glmer()公式还是应该传递一些指针和引用？我在示例中看到X(bootMer()的第一个参数是*lmer()对象.我想编写*Mixed5但是它呈现错误.

非常感谢.

我的代码是:

library(lme4)
library(boot)

(mixed5 <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
                 + (1 | PatientID) + (0 + Trt | PatientID)
                 , family=binomial(logit), MixedModelData4))


FUN <- function(formula) {
  fit <- glmer(DV ~ (Demo1 +Demo2 +Demo3 +Demo4 +Trt)^2 
               + (1 | PatientID) + (0 + Trt | PatientID)
               , family=binomial(logit), MixedModelData4)
  return(coef(fit))
}

result <- bootMer(mixed5, FUN, nsim = 3, seed = NULL, use.u = FALSE,
        type = c("parametric"),
        verbose = T, .progress = "none", PBargs = …