我创建了一个类似下面的脚本来做我称之为"加权"回归的事情:
library(plyr)
set.seed(100)
temp.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)
temp.df.expand <- ddply(temp.df,
c("uid"),
function(df) {
data.frame(bp=rep(df[,"bp"],df[,"weight"]),
age=rep(df[,"age"],df[,"weight"]),
stringsAsFactors=FALSE)})
temp.df.lm <- lm(bp~age,data=temp.df,weights=weight)
temp.df.expand.lm <- lm(bp~age,data=temp.df.expand)
Run Code Online (Sandbox Code Playgroud)
你可以看到,在temp.df,每一行都有它的重量,我的意思是,共有1178样本,但对于相同的行为bp及age,它们是合并成1行的代表weight列.
我使用了weight函数中的参数lm,然后用另一个数据帧交叉检查结果,数据temp.df帧是"扩展"的.但我发现lm2个数据帧的输出不同.
我是否误解了weight函数中的参数lm,并且任何人都可以告诉我如何正确运行回归(即不手动扩展数据帧)以获得类似的数据集temp.df吗?谢谢.
我的首要问题是,R如何在WLS案例中计算R ^ 2?它不仅仅对观测值进行加权,然后计算R ^ 2.为了解决这个问题,我正在浏览源代码,直到我在代码中遇到这个问题lm.wfit:
z <- .Call(C_Cdqrls, x *wts, y*wts, tol)
Run Code Online (Sandbox Code Playgroud)
在这做什么?有谁知道如何访问代码以获取详细信息?即,返回的是z什么?如何C_Cdqrls,x*wts,y*wts,tol被使用?
到目前为止我所理解的(我不确定它是否正确),这.Call意味着R在C中执行此代码.但是,如果可能的话,我想看看如何在C中完成此操作.
谢谢!
我对两个值进行了数十亿次测量,x并且y. 这太大而无法对原始数据进行操作,因此我将它们表示为频率表。对于每个唯一的x值和y值组合,我都有一行,还有一个变量freq显示有多少数据点具有该值的组合。
如果我想估计x和y之间的关系,我可以这样做:lm(y ~ x, data=df, weights=df$freq)。我已经对此进行了测试,它给出了准确的参数估计,但给出了错误的t值。它仍然将每一行视为一个观察,因此自由度比它们应该的要小得多。
注意:这个问题显示了如何重新创建原始数据,但我的原始数据非常大,这就是我首先使用频率表的原因。
例子
# This dataset has a negative correlation between x and y:
library(dplyr)
raw_data<-data.frame(
x=rep(c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4), 100),
y=rep(c(5,5,5,5,1,4,4,4,4,1,3,3,3,3,7,2,2,2,2,8), 100)
)
lm_raw<-lm(x ~ y, data=raw_data)
summary(lm_raw)[c("coefficients", "df")]
# Let's say instead I have a have a summary dataset that has the frequency for each x-y pair:
freq_data <- raw_data %>% group_by(x,y) %>% summarise(freq=n())
# …Run Code Online (Sandbox Code Playgroud)