r - 使用带有 lm 的样本权重时更正 dfs
问题描述
我试图弄清楚权重lm
实际上是如何工作的,我看到了这个 7.5 年前的问题,它提供了一些关于权重如何工作的见解。该问题的数据部分复制并在下面扩展。
我在 Cross Validated 上发布了这个相关的问题。
library(plyr)
set.seed(100)
df <- data.frame(uid=1:200,
bp=sample(x=c(100:200),size=200,replace=TRUE),
age=sample(x=c(30:65),size=200,replace=TRUE),
weight=sample(c(1:10),size=200,replace=TRUE),
stringsAsFactors=FALSE)
set.seed(100)
df.double_weights <- data.frame(uid=1:200,
bp=sample(x=c(100:200),size=200,replace=TRUE),
age=sample(x=c(30:65),size=200,replace=TRUE),
weight=2*df$weight,
stringsAsFactors=FALSE)
df.expand <- ddply(df,
c("uid"),
function(df) {
data.frame(bp=rep(df[,"bp"],df[,"weight"]),
age=rep(df[,"age"],df[,"weight"]),
stringsAsFactors=FALSE)})
df.lm <- lm(bp~age,data=df,weights=weight)
df.double_weights.lm <- lm(bp~age,data=df.double_weights,weights=weight)
df.expand.lm <- lm(bp~age,data=df.expand)
summary(df.lm)
summary(df.double_weights.lm)
summary(df.expand.lm)
这三个 data.frames 包含完全相同的数据。然而;
其中df
有 200 个观察值,加权加起来为 1178 sum(df.$weight) == 1178
,.
在df.double_weights
中,权重只是翻了一番sum(df.double_weights$weight) == 2356
。
在df.expand
中,有 1178 个未加权观测值,而不是 200 个加权观测值。
两者的系数summary(df.lm)
相同summary(df.double_weights.lm)
,重要性也相同(这意味着,如果权重工作正常,权重的绝对大小无关紧要)。编辑:但是,绝对大小似乎确实很重要,请参阅底部结果。
但是,对于summary(df.lm)
和summary(df.expand.lm)
,系数相同,但显着性不同。
summary(df.lm)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 165.6545 10.3850 15.951 <2e-16 ***
age -0.2852 0.2132 -1.338 0.183
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 98.84 on 198 degrees of freedom
Multiple R-squared: 0.008956, Adjusted R-squared: 0.003951
F-statistic: 1.789 on 1 and 198 DF, p-value: 0.1825
summary(df.expand.lm)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 165.65446 4.26123 38.88 < 2e-16 ***
age -0.28524 0.08749 -3.26 0.00115 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 28.68 on 1176 degrees of freedom
Multiple R-squared: 0.008956, Adjusted R-squared: 0.008114
F-statistic: 10.63 on 1 and 1176 DF, p-value: 0.001146
根据@IRTFM 的说法,自由度没有被正确添加,提供了以下代码来修复它:
df.lm.aov <- anova(df.lm)
df.lm.aov$Df[length(df.lm.aov$Df)] <-
sum(df.lm$weights)-
sum(df.lm.aov$Df[-length(df.lm.aov$Df)] ) -1
df.lm.aov$`Mean Sq` <- df.lm.aov$`Sum Sq`/df.lm.aov$Df
df.lm.aov$`F value`[1] <- df.lm.aov$`Mean Sq`[1]/
df.lm.aov$`Mean Sq`[2]
df.lm.aov$`Pr(>F)`[1] <- pf(df.lm.aov$`F value`[1], 1,
df.lm.aov$Df, lower.tail=FALSE)[2]
df.lm.aov
Analysis of Variance Table
Response: bp
Df Sum Sq Mean Sq F value Pr(>F)
age 1 8741 8740.5 10.628 0.001146 **
Residuals 1176 967146 822.4
现在,将近 8 年后,显然这个问题仍然存在(这是否意味着几乎所有使用加权变量与lm
from结合的研究的R
显着性值都太低了?)更实际地,我遇到的问题是我几乎不明白 IRTFM 是什么正在做什么,或者它与多元回归分析(甚至是在后台使用lm
的其他函数)有何关系?
问题:有没有更通用的方法来解决这个问题,可以应用于多元回归?
编辑:
如果我们在 上运行 IRTFM 的解决方案df.double_weights.lm
,我们会得到不同的结果,因此显然权重的绝对大小确实很重要。
Analysis of Variance Table
Response: bp
Df Sum Sq Mean Sq F value Pr(>F)
age 1 17481 17481.0 21.274 4.194e-06 ***
Residuals 2354 1934293 821.7
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
解决方案
如果我正确理解您的问题,那么您在权重列中的内容通常称为“频率权重”。它们用于通过指示您对每个协变量组合有多少观察值来节省数据集中的空间。
要使用“聚合”数据集估计模型并获得正确的标准误差,您需要做的就是更正模型中的自由度数lm
。
正确的自由度数是观察总数减去模型中的参数数。这可以通过获取weights
变量的总和或通过查看“完整”数据中的观察总数并减去估计的参数数量(即系数)来计算。
这是一个更简单的例子,我认为这更清楚了这一点:
library(dplyr)
library(modelsummary)
set.seed(1024)
# individual (true) dataset
x <- round(rnorm(1e5))
y <- round(x + x^2 + rnorm(1e5))
ind <- data.frame(x, y)
# aggregated dataset
agg <- ind %>%
group_by(x, y) %>%
summarize(freq = n())
models <- list(
"True" = lm(y ~ x, data = ind),
"Aggregated" = lm(y ~ x, data = agg),
"Aggregated & W" = lm(y ~ x, data = agg, weights=freq),
"Aggregated & W & DF" = lm(y ~ x, data = agg, weights=freq)
)
现在我们要更正列表中最后一个模型的自由度数。我们通过对我们的freq
列求和来做到这一点。我们也可以使用nrow(ind)
,因为它们是相同的:
# correct degrees of freedom
models[[4]]$df.residual <- sum(agg$freq) - length(coef(models[[4]]))
最后,我们总结了使用该modelsummary
包的所有 5 个模型。请注意,第一个模型和最后一个模型完全相同,即使第一个模型是使用完整的单个数据集估计的,而最后一个模型是使用聚合数据估计的:
modelsummary(models, fmt=5)
真的 | 聚合的 | 聚合 & W | 聚合 & W & DF | |
---|---|---|---|---|
(截距) | 1.08446 | 5.51391 | 1.08446 | 1.08446 |
(0.00580) | (0.71710) | (0.22402) | (0.00580) | |
X | 1.00898 | 0.91001 | 1.00898 | 1.00898 |
(0.00558) | (0.30240) | (0.21564) | (0.00558) | |
编号Obs。 | 1e+05 | 69 | 69 | 69 |
R2 | 0.246 | 0.119 | 0.246 | 0.246 |
R2 调整。 | 0.246 | 0.106 | 0.235 | 0.999 |
AIC | 405058.1 | 446.0 | 474.1 | 474.1 |
BIC | 405086.7 | 452.7 | 480.8 | 480.8 |
日志.Lik。 | -202526.074 | -219.977 | -234.046 | -234.046 |
F | 32676.664 | 9.056 | 21.894 | 32676.664 |
推荐阅读
- reactjs - React Typescript 不断收到无效的钩子调用,不知道如何解决?
- java - 如何将过滤器集合应用于 Java 流?
- ruby-on-rails - RAILS - 如何编写一个查询来获取彼此在 1 分钟内创建的记录
- python - python中用户定义似然的贝叶斯分析
- java - 捕获到意外异常:运行 AsyncInstallTest 时出现 com.ptc.pfc.Implementation.pfcExceptions$XToolkitGeneralError
- hadoop - 启用 tez 后 Hive 查询总是进行故障转移
- angular - 这是对 Observables 进行排队的合理方式吗?
- flutter - 升级到空安全但未定义接收“地图”
- python - 我有一个 Selenium 脚本,可以在 Chrome 中打开一个 URL 并将网页保存为 PDF。我如何重写它以在 Safari 中做同样的事情?
- mongodb - Hazelcast MapStore 的 store() 方法在 map.put() 操作中没有被调用,在客户端/服务器模式下,但它在嵌入式模式下完美运行