0

我根据三个月的时间段对一些数据框进行了子集化，并命名为 jfm（1 月至 3 月）、fma（2 月至 4 月）、mam（3 月至 5 月）……直到 ond（10 月至 12 月）。我希望使用几个变量作为回归量对所有这些数据进行类似的分析。下面我将展示如何使用其中一种污染物作为回归量对两个子集数据帧之一进行分析。我有兴趣对分别输入模型的所有污染物（pm10median、pm25median、o3median 和 so2median）进行分析。如何对所有数据框进行此分析？

library(gamair) 
library(mgcv)
data(chicago) 
chicago$date<-seq(from=as.Date("1987-01-01"), to=as.Date("2000-12-31"),length=5114)


chicago$month<-as.numeric(format(chicago$date,"%m")) ## create month
jfm <- subset(chicago, month %in% c(1:3) )      ## subset data for January to March
fma <- subset(chicago, month %in% c(2:4) )  ## February to April
mam <- subset(chicago, month %in% c(3:5) )  ## March to may


jfm$trend<-seq(dim(jfm)[1])   ## cretae a trend for specific df based on dimension of the df
fma$trend<-seq(dim(fma)[1])   ## trend for df fma


## Regress each pollutant separately on death for the first subset

model1<-gam(death ~  pm10median + s(trend,k=21)+ s(tmpd,k=6) ,family=quasipoisson,na.action=na.omit,data=jfm) 

model2<-gam(death ~  pm10median + s(trend,k=21)+ s(tmpd,k=6) ,family=quasipoisson,na.action=na.omit,data=fma)
4

1 回答 1

0 
# create a function that defines the exact regression
# you want to run on all three-month data sets
fun <- 
    function( y , x ){

        # store each of the regression outputs into an object
        a <- gam(
            death ~  pm10median + s(trend,k=21)+ s(tmpd,k=6) ,
            family = quasipoisson , 
            na.action = na.omit ,
            data = x[ x$month %in% y , ]
        ) 
        b <- gam(
            death ~  pm25median + s(trend,k=21)+ s(tmpd,k=6) ,
            family = quasipoisson , 
            na.action = na.omit ,
            data = x[ x$month %in% y , ]
        ) 

        # return each of the regressions as a list
        list( a , b )
    }

# define which three-month groups you want to run it on
months <- cbind( 1:10 , 2:11 , 3:12 )

# now just run the function for each row in `months`
results <- apply( months , 1 , fun , x = chicago )

# look at the whole thing
results

# extract jfm, for example
jfm <- results[[1]]

# extract fma (and print it to the screen as well)
( fma <- results[[2]] )
于 2013-02-05T09:36:27.973 回答