R：引导多重回归答案

【问题标题】：R: Bootstrap Multiple RegressionR：引导多重回归
【发布时间】：2019-03-02 18:59:56
【问题描述】：

我正在尝试进行多元回归分析，以找出水质对特定位置（又名 Guzzler）的浮游生物丰度的影响。我能够让我的模型运行和汇总，但是数据是非参数的，因此典型的汇总是不可靠的。这主要是因为样本量很小，因为它是在几周的过程中完成的，每周一个样本。

当时我在想这个的非参数版本可能是一个引导程序。我以前在其他数据上运行过引导程序，但从来没有使用多元回归模型。我似乎找不到关于如何解决这个问题的代码，所以从我过去如何执行引导程序开始。我很好奇我需要编辑什么才能让这个引导程序运行。

这是 dput(head(Guzzler1) 的输出：

    structure(list(Abundance = c(98L, 43L, 65L, 55L, 54L), Phospates = c(2L, 
2L, 2L, 2L, 2L), Nitrates = c(0, 0.3, 0, 0.15, 0), pH = c(7.5, 
8, 7.5, 7, 7)), .Names = c("Abundance", "Phospates", "Nitrates", 
"pH"), row.names = c(NA, 5L), class = "data.frame")

这是我的模型和摘要：

    Guzzler1model<-lm(Abundance ~ Phospates + Nitrates + pH, data=Guzzler1)
> summary(Guzzler1model)

Call:
lm(formula = Abundance ~ Phospates + Nitrates + pH, data = Guzzler1)

Residuals:
     1      2      3      4      5 
 20.75  -4.25 -12.25   8.50 -12.75 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   -80.25     209.62  -0.383    0.739
Phospates         NA         NA      NA       NA
Nitrates     -135.00      90.02  -1.500    0.272
pH             21.00      28.87   0.727    0.543

Residual standard error: 20.41 on 2 degrees of freedom
Multiple R-squared:  0.5302,    Adjusted R-squared:  0.06032 
F-statistic: 1.128 on 2 and 2 DF,  p-value: 0.4698

***请注意：我相信磷酸盐具有 NA，因为在这个特定位置每个值都等于 2。

这是我最初执行引导程序但不确定要更改什么的方式：

n=length(Guzzler1Abundance)
B = 1000
results = numeric(B)
for(b in 1:B){
  i = sample(x = 1:n, size=n, replace=TRUE)
  bootSample = Guzzler1Abundance[i]
  thetahat= mean(bootSample)
  results[b] = thetahat
}

提前非常感谢您！

【问题讨论】：

我认为您需要包含更大的数据集以支持引导演示。 5 个四分量模型的案例很可能会出错。
@42- 不幸的是，我无法包含更多数据，因为这是针对某些学生的项目。您对我运行的回归模型的非参数版本有什么建议吗？
这是一个上帝的起点：socialsciences.mcmaster.ca/jfox/Books/Companion/appendices/…

标签： r model linear-regression statistics-bootstrap

【解决方案1】：

我不太清楚你所说的非参数数据是什么意思，但我理解，你想从你的数据中提取自举样本并用它执行线性回归。

一种可能的方法是

Guzzler1 <-  structure(list(Abundance = c(98L, 43L, 65L, 55L, 54L), Phospates = c(2L, 
                      2L, 2L, 2L, 2L), Nitrates = c(0, 0.3, 0, 0.15, 0), pH = c(7.5, 
                      8, 7.5, 7, 7)), .Names = c("Abundance", "Phospates", "Nitrates", 
                      "pH"), row.names = c(NA, 5L), class = "data.frame")
lines <- nrow(Guzzler1)
replicate(5, lm(Abundance ~ Phospates + Nitrates + pH, 
                         data=Guzzler1[sample(lines, replace = TRUE),])$coefficients)

这将报告来自像这样的 5 个线性回归的系数

> replicate(5, lm(Abundance ~ Phospates + Nitrates + pH, 
+                          data=Guzzler1[sample(lines, replace = TRUE),])$coefficients)
                 [,1]       [,2]      [,3]       [,4] [,5]
(Intercept)  65.00000 145.000000 -408.0000 145.000000 -100
Phospates          NA         NA        NA         NA   NA
Nitrates    -73.33333   6.666667 -256.6667   6.666667 -110
pH                 NA -13.000000   66.0000 -13.000000   22

通过更改我的replicate 调用的第一个参数，可以任意选择更高的重复次数。正如@IRTFM 预测和解释的那样，许多NA 值是由于数据稀缺造成的。随着对更多数据的采样，这种情况会有所改善。

让我们抽取 5000 个 bootstrap 样本并研究 Nitrate 系数的分布：

reps <- replicate(5000, lm(Abundance ~ Phospates + Nitrates + pH, 
                           data=Guzzler1[sample(lines, replace = TRUE),])$coefficients)
plot(table(reps["Nitrates",]))
plot(ecdf(reps["Nitrates",]))
quantile(reps["Nitrates",], c(.025, .25, .5, .75, .975), na.rm = TRUE)

磷酸盐可以编辑到这个，一有变异数据：

boxplot(reps["(Intercept)",], reps["Nitrates",], reps["pH",]
        , names = c("Intercept", "Nitrates", "pH"), ylab="bootstrapped coefficients")
abline(h=0, col="firebrick", lty=3)

【讨论】：