2

我有一个包含多个变量的大型数据集,其中一个是状态变量,每个状态编码为 1-50。我想对数据集的剩余 27 个变量(总共 55 个变量)进行 28 个变量的回归,并且针对每个州进行回归。

换句话说,对变量 1 在 covariate1、covariate2、...、covariate27 上运行回归,以观察 state==1。然后,我想对状态 2-50 的变量 1 重复此操作,并对变量 2、变量 3、...、变量 28 重复整个过程。

我想我已经编写了正确的 R 代码来执行此操作,但接下来我想做的是提取系数,最好是提取到系数矩阵中。有人可以帮我吗?这是我到目前为止编写的代码:

for (num in 1:50) {

    #PUF is the data set I'm using

    #Subset the data by states
    PUFnum <- subset(PUF, state==num)

    #Attach data set with state specific data
    attach(PUFnum)

    #Run our prediction regression
    #the variables class1 through e19700 are the 27 covariates I want to use
    regression <- lapply(PUFnum,  function(z) lm(z ~ class1+class2+class3+class4+class5+class6+class7+
                                                     xtot+e00200+e00300+e00600+e00900+e01000+p04470+e04800+
                                                     e09600+e07180+e07220+e07260+e06500+e10300+
                                                     e59720+e11900+e18425+e18450+e18500+e19700))

    Beta <- lapply(regression, function(d) d<- coef(regression$d))


    detach(PUFnum)
}
4

2 回答 2

5

这是另一个经典Split-Apply-Combine问题的例子,可以使用plyr@hadley 的包来解决。在你的问题中,你想

  1. 按状态拆分数据帧
  2. 对每个子集应用回归
  3. 将系数组合到数据框中。

我将用图书馆中Cars93可用的数据集来说明它。MASS我们有兴趣弄清楚国家之间的关系以及horsepower国家之间的关系。enginesizeorigin

# LOAD LIBRARIES
require(MASS); require(plyr)

# SPLIT-APPLY-COMBINE
regressions <- dlply(Cars93, .(Origin), lm, formula = Horsepower ~ EngineSize)
coefs <- ldply(regressions, coef)

   Origin (Intercept) EngineSize
1     USA    33.13666   37.29919
2 non-USA    15.68747   55.39211

编辑。对于您的示例,替换PUFCars93, stateforOriginfmfor the formula

于 2012-01-26T20:35:17.047 回答
1

我已经稍微清理了你的代码:

fm <- z ~ class1+class2+class3+class4+class5+class6+class7+
          xtot+e00200+e00300+e00600+e00900+e01000+p04470+e04800+
          e09600+e07180+e07220+e07260+e06500+e10300+
          e59720+e11900+e18425+e18450+e18500+e19700

PUFsplit <- split(PUF, PUF$state)
mod <- lapply(PUFsplit, function(z) lm(fm, data=z))

Beta <- sapply(mod, coef)

如果您愿意,您甚至可以将所有内容放在一行中:

Beta <- sapply(lapply(split(PUF, PUF$state), function(z) lm(fm, data=z)), coef)
于 2012-01-26T06:47:56.077 回答