使用带有 lm 的样本权重时更正 dfs
Correcting dfs when using sample weights with lm
我试图弄清楚 lm 中的权重实际上是如何工作的，我看到了 this 7,5 year old question，它提供了一些关于权重如何工作的见解。该问题的数据部分复制并在下面扩展。
我在 Cross Validated 上发布了 this related question。
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 年过去了，显然这个问题仍然存在（这是否意味着几乎所有使用加权变量与 R 中的 lm 相结合的研究都具有太低的显着性值？ ) 更实际地说，我遇到的问题是我几乎不了解 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
如果我正确理解你的问题，你在权重列中的内容通常称为“频率权重”。它们用于通过指示每个协变量组合有多少观察值来将 space 保存在数据集中。
要使用“聚合”数据集估计模型并获得正确的标准误差，您需要做的就是更正 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)
True
Aggregated
Aggregated & W
Aggregated & W & DF
(Intercept)
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)
Num.Obs.
1e+05
69
69
69
R2
0.246
0.119
0.246
0.246
R2 Adj.
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
Log.Lik.
-202526.074
-219.977
-234.046
-234.046
F
32676.664
9.056
21.894
32676.664
我试图弄清楚 lm 中的权重实际上是如何工作的，我看到了 this 7,5 year old question，它提供了一些关于权重如何工作的见解。该问题的数据部分复制并在下面扩展。
我在 Cross Validated 上发布了 this related question。
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 年过去了，显然这个问题仍然存在（这是否意味着几乎所有使用加权变量与 R 中的 lm 相结合的研究都具有太低的显着性值？ ) 更实际地说，我遇到的问题是我几乎不了解 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
如果我正确理解你的问题，你在权重列中的内容通常称为“频率权重”。它们用于通过指示每个协变量组合有多少观察值来将 space 保存在数据集中。
要使用“聚合”数据集估计模型并获得正确的标准误差，您需要做的就是更正 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)
| True | Aggregated | Aggregated & W | Aggregated & W & DF | |
|---|---|---|---|---|
| (Intercept) | 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) | |
| Num.Obs. | 1e+05 | 69 | 69 | 69 |
| R2 | 0.246 | 0.119 | 0.246 | 0.246 |
| R2 Adj. | 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 |
| Log.Lik. | -202526.074 | -219.977 | -234.046 | -234.046 |
| F | 32676.664 | 9.056 | 21.894 | 32676.664 |