使用sapply函数改写for循环并绘制迟滞温度与污染物效应图

For循环应该是我们在R语言使用得最普遍的循环了，优势就是简单、易上手，缺点就是慢，特别对于跑数据量比较大的数据。Apply家族函数使用C++来编写，运行得非常快，非常适合代替for循环。今天介绍一下sapply函数改写for循环并绘制迟滞温度与污染物效应图。
我们先导入数据和R包

library(mgcv)
library(splines)
library(tsModel)
bc<-read.csv("E:/r/test/chicago.csv",sep=',',header=TRUE)

这是美国芝加哥1987年至 2000年大气污染与死亡数据（公众号回复：芝加哥2，可以获得数据），我们先来看看数据的构成，death：死亡人数 (per day)，pm10：大气污染物pm10的中位数值，pm25median，o3median：二氧化硫的中位数值，time：天数，这里就是我们的时间，tmpd：华氏温度，date：日期
假设我们研究温度的迟滞和死亡率关系，研究2周内的迟滞关系，需要写个循环
先生成一个向量和一个空值，等下跑循环用

maxlag <- 0:13
b<-NULL
for (i in maxlag) {
  fit <- glm(death~pm10 + Lag(temp,seq(0, i)), data = bc, family = poisson)
  summ <- summary(fit)
  summ.coef <- summ$coefficients["pm10",2]
  c(coef(fit)["pm10"], summ.coef)
  a<-rbind(coef(fit)["pm10"],summ.coef)
  b<-cbind(b,a)
}
models<-b

这样就把2周的温度迟滞对死亡率的影响的数据跑出来了，我们画个图

rng <- range(models[1, ] - 1.96*models[2,], models[1, ] + 1.96*models[2,], 0)
par(mar = c(4, 5, 1, 1))
plot(maxlag, models[1, ], type = "b",pch = 20, ylim = rng, xlab = "Maximum temperature lag",
     ylab = expression(hat(beta)*" for "* PM[10]*" at lag 1"))
lines(maxlag, models[1, ] + 1.96*models[2, ], lty = 2)
lines(maxlag, models[1, ] - 1.96*models[2,], lty = 2)
abline(h = 0, lty = 3)

这样，温度迟滞对死亡率关系的影响曲线就出来了。下面我们通过sapply函数来改写这个for循环
先要写一个程序来解释一下这个for循环

function(mlag) {
  fit <- glm(death~pm10 + Lag(temp,seq(0, mlag)), data = bc, family = poisson)
  summ <- summary(fit)
  summ.coef <- summ$coefficients["pm10",2]
  c(coef(fit)["pm10"], summ.coef)
  }

写好以后使用sapply函数把它包起来就可以了

models <- sapply(maxlag, function(mlag) {
  fit <- glm(death~pm10 + Lag(temp,seq(0, mlag)), data = bc, family = poisson)
  summ <- summary(fit)
  summ.coef <- summ$coefficients["pm10",2]
  c(coef(fit)["pm10"], summ.coef)
  })

两个循环生成的models函数一模一样

最后绘图也是一样的，我这个数据小，看不出来差别，大型数据跑起来明显不一样的。