使用 by() 从模型中提取预测值

Question

我正在将指数模型拟合到 208 个泉水区的人口数据中，以 5 年为间隔回溯计算 1975-2015 年的人口，即 seq(1975,2015,5) .这是我数据集中的前 5 个弹簧，以及我用来拟合模型并绘制它的代码（我想有数字）：

springsheds <- 
structure(list(spring = c("alexander", "alexander", "alexander", "alexander", 
"blue hole", "blue hole", "blue hole", "blue hole", "cedar head", "cedar head", 
"cedar head", "cedar head", "charles", "charles", "charles", "charles", 
"columbia", "columbia", "columbia", "columbia"), year = c(2000L, 2005L, 2010L, 
2015L, 2000L, 2005L, 2010L, 2015L, 2000L, 2005L, 2010L, 2015L, 2000L, 2005L, 
2010L, 2015L, 2000L, 2005L, 2010L, 2015L), pop = c(527L, 620L, 732L, 867L, 
3071L, 3356L, 3669L, 4007L, 3038L, 3320L, 3630L, 3965L, 1311L, 1446L, 1592L, 
1747L, 7550L, 8130L, 8706L, 9332L)), .Names = c("spring", "year", "pop"), 
class = "data.frame", row.names = c(NA, -20L))
models.spsh <- by(springsheds, springsheds$spring, function(x) {
        fm <- lm(log(pop) ~ year, data = x)
        timevalues <- seq(1970, 2020, 10)
        predict <- exp(predict(fm,list(year=timevalues)))
        plot(pop ~ year, x, main = spring[1], xlim = c(1970, 2020), ylim=c(0,15000))
        lines(timevalues, predict,lwd=1, col = "blue", xlab = "Year", ylab = "Population")
})

我也可以使用 by（） 来提取每个弹簧的预测值吗？我目前的解决方法是为每个弹簧单独创建一个对象，并将预测值迭代添加到一个对象中：

fm <- lm(log(pop) ~ year, data = alex)
timevalues <- seq(1975,2015,5)
alex <- exp(predict(fm,list(year=timevalues)))
old<-cbind(timevalues,alex)
fm <- lm(log(pop) ~ year, data = blue)
blue <- exp(predict(fm,list(year=timevalues)))
old<-cbind(old,blue)

这似乎效率很低，我假设有一种更优雅的方法可以做到这一点，有没有一种方法可以添加到我的初始代码中以提取预测的总体值？

Answer 1

您可以split数据，然后对每个所需的输出使用 lapply：

#Split the data grouped by spring
sdata <- split(springsheds, springsheds$spring)
#Fit the models
fit.spsh <- lapply(sdata, function(x) {
  lm(log(pop) ~ year, data = x)
})
#Get the predicted values
timevalues <- seq(1970, 2020, 10)
predictList <- lapply(fit.spsh, function(m) exp(predict(m,list(year=timevalues))))
#Generate plots
lapply(names(sdata), function(n) {
  plot(pop ~ year,sdata[[n]] , main = n, xlim = c(1970, 2020), ylim=c(0,15000))
  lines(timevalues, predictList[[n]],lwd=1, col = "blue", xlab = "Year", ylab = "Population")
})
 #Combine the predict values
 do.call(cbind,predictList)
 #alexander blue hole cedar head   charles  columbia  
 #1  194.3679  1803.470   1783.068  738.9545  4955.633
 #2  270.8778  2153.663   2129.682  894.9253  5705.076
 #3  377.5048  2571.856   2543.676 1083.8167  6567.857
 #4  526.1037  3071.253   3038.146 1312.5774  7561.118
 #5  733.1965  3667.621   3628.738 1589.6225  8704.590
 #6 1021.8081  4379.790   4334.137 1925.1434 10020.989