14

「加重」回帰と呼ばれるものを実行するために、以下のようなスクリプトを作成しました。

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)

temp.dfでは、各行に重みがあることがわかります。つまり、合計 1178 のサンプルがありますが、同じbpとの行でageは、それらは 1 行にマージされ、weight列に表示されます。

weight関数でパラメーターを使用した後、データフレームが「展開」されlmている別のデータフレームと結果をクロスチェックします。temp.dfしかしlm、2つのデータフレームの出力が異なることがわかりました。

weight関数内のパラメーターを誤って解釈しましたか?lmのように提示されたデータセットに対して回帰を適切に (つまり、データフレームを手動で拡張せずに) 実行する方法を誰かに教えてもらえますtemp.dfか? ありがとう。

4

1 に答える 1

16

ここでの問題は、適切な Df 統計と平均平方和統計を取得するために、自由度が適切に加算されていないことです。これにより、問題が修正されます。

temp.df.lm.aov <- anova(temp.df.lm)
temp.df.lm.aov$Df[length(temp.df.lm.aov$Df)] <- 
        sum(temp.df.lm$weights)-   
        sum(temp.df.lm.aov$Df[-length(temp.df.lm.aov$Df)]  ) -1
temp.df.lm.aov$`Mean Sq` <- temp.df.lm.aov$`Sum Sq`/temp.df.lm.aov$Df
temp.df.lm.aov$`F value`[1] <- temp.df.lm.aov$`Mean Sq`[1]/
                                        temp.df.lm.aov$`Mean Sq`[2]
temp.df.lm.aov$`Pr(>F)`[1] <- pf(temp.df.lm.aov$`F value`[1], 1, 
                                      temp.df.lm.aov$Df, lower.tail=FALSE)[2]
temp.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

と比べて：

> anova(temp.df.expand.lm)
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                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

これが R-help で頻繁に出てこないことに少し驚いています。それか、私の検索戦略の開発力が年をとって弱まっています。

于 2012-04-22T18:47:57.077 に答える