R中回归的ddply

Question

我有一个数据框架，其中包含了100个城市的“输出”、平均温度、湿度和时间(作为24个因素，而不是连续的)数据(由代码提供)。我想应用一个回归公式，根据温度、湿度和时间数据来预测每个城市的产量。我希望得到100个不同的回归模型。我使用了ddply函数，并在this thread的帮助下编写了下面的代码行。

df = ddply(data, "city", function(x) coefficients(lm(output~temperature+humidity, data=x)))

此代码适用于数值数据、温度和湿度。但是，当我添加时区因素数据(也就是23个因子变量)时，我会得到一个错误：

df = ddply(data, "city", function(x) coefficients(lm(output~temperature+humidity+time, data=x)))

错误:对比只能适用于2级或2级以上的因素

有人知道这是为什么吗？下面是我的数据框架的一个示例块：

city    temperature   humidity   time   output
 11        51            34        01     201
 11        43            30        02     232
 11        55            50        03     253  
 11        64            54        10     280  
 22        21            52        11     321  
 22        43            65        04     201  
 22        51            66        09     211  
 22        51            78        16     199  
 05        45            70        01     202  
 05        51            54        10     213 

所以我想要三个城市的模型，基于温度，湿度和时间因素。

Answer 1

通过使用ddply，可以将lm应用于数据帧的子集，其中每个子集对应于某个城市。在完整的数据集中，有些城市似乎只有一个记录。对于这种情况，统计分析显然是毫无意义的，但是lm会返回一些答案，但是如果模型中有一个因子变量，它会抛出一个错误。

作为解决办法，您可以检查匿名函数中的行数：

ddply(d,'city',function(x) if (nrow(x)==1) return() else coefficients(lm(output~temperature+humidity+time, data=x)))

其中d是您的示例集的稍微修改的版本，在该版本中，我在最后一行更改了城市的id，以确保某些城市只有一条记录：

d <- structure(list(city = c(11, 11, 11, 11, 22, 22, 22, 22, 5, 7), temperature = c(51L, 43L, 55L, 64L, 21L, 43L, 51L, 51L, 45L,     51L), humidity = c(34L, 30L, 50L, 54L, 52L, 65L, 66L, 78L,     70L, 54L), time = structure(c(1L, 2L, 3L, 6L, 7L, 4L, 5L,     8L, 1L, 6L), .Label = c("1", "2", "3", "4", "9", "10", "11",     "16"), class = "factor"), output = c(201L, 232L, 253L, 280L,     321L, 201L, 211L, 199L, 202L, 213L)), .Names = c("city", "temperature", "humidity", "time", "output"), row.names = c(NA, -10L), class = "data.frame")

您也可以使用这个R基代码而不是ddply

L <- split(d,d$city)

L2 <- lapply(L,function(x) {
    if (nrow(x)==1) 
        return() 
    else 
        coefficients(lm(output~temperature+humidity+time, data=x))
})

M <- do.call(rbind,L2)
df <- as.data.frame(M)

这段代码更冗长，但在出现问题行为时，检查和分析它要容易得多。