在 R 中的 for 循环中运行时间序列滞后

Question

我遇到了计算和运行降水变异系数滞后的 for 循环问题。我不太确定如何概括这个问题，所以我添加了到目前为止我已经采取的所有步骤。

我的主要数据集“d”如下所示：

      row.names timestamp  station      year   month   ndvi  landcover    altitude precipitation
1         1         1      A            2000   jan   0.4138  Mixed forest     2143          16.0
2      1769         2      A            2000   feb   0.4396  Mixed forest     2143           4.0

我想找出降水变异系数滞后 0:10 对每个站点一年的最大 ndvi 的影响。 基本上我的代码如下所示：

r <- aggr(d,c("station","landcover","year"), c("altitude=mean(altitude)","max.ndvi=NA","max.month=NA","max.timestamp=NA","max.precipitation=NA", "cv=NA"))


head(r)
    station    landcover year altitude max.ndvi max.month max.timestamp max.precipitation cv
1     A      Mixed forest 2000     2143       NA        NA            NA          NA     NA
2     A      Mixed forest 2001     2143       NA        NA            NA          NA     NA

for(i in 1:nrow(r)) {
  tmp <- d[d$station==r$station[i] & d$year==r$year[i],]
  idx <- which.max(tmp$ndvi);
  r$max.month[i] <- as.character(tmp$month[idx]);   
  r$max.ndvi[i] <- tmp$ndvi[idx];
  r$max.timestamp[i] <- tmp$timestamp[idx];
  r$max.precipitation[i] <- tmp$precipitation[idx];
  r$cv[i] <- sd(tmp$precipitation, na.rm = TRUE)/mean(tmp$precipitation, na.rm = TRUE) 
}

for(lag in 0:10) {
  cat("\n\n***** lag =",lag,"*****\n\n");
  for(i in 1:nrow(r)) {
    timestamp <- r$max.timestamp[i]-lag;
    if(timestamp>0){
    r$cv[i] <- r$cv[d$station==r$station[i] & d$timestamp==timestamp];
    }
  }
  r <- na.omit(r)
  print(summary(aov(max.ndvi~cv, data=r)));
  for(lu in sort(unique(as.character(r$landcover)))) {
  cat("\n----------------- Analysis for LU =",lu,"\n\n");
  print(summary(aov(max.ndvi~cv,data=r[r$landcover==lu,])));
  }
}

我遇到的问题是最后一部分为每个 max.ndvi 值分配/循环滞后。我想要所有行的每个滞后的摘要以及每种土地覆盖类型的摘要。

我尝试了各种不同的组合，但总是出错。对于上面的代码，我得到了这个错误：

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  0 (non-NA) cases 

谁能给点建议？

非常感谢。

Answer 1

我认为您的一些 cv lag 有所有缺失值（在土地覆盖的详细程度）。此代码要求 cv.lags（在土地覆盖类中）至少有 3 个观测值。

#Create fake dataset in your same mold
set.seed(1)
d <- data.frame(row.names=1:40,timestamp=rep(1:10,4),station=c(rep("A",10),rep("B",10),rep("C",10),rep("D",10)),
                year=rep(c(2000,2000,2000,2001,2001,2001,2001,2002,2002,2002),4),
                month=rep(c("jan","feb","mar","april","may","jun","jul","aug","sep","oct"),4),ndvi=runif(40,min=0.3,max=0.6),
                landcover=c(rep("Mixed Forest",10),rep("Sand",10),rep("other1",10),rep("other2",10)),altitude=runif(40,min=1500,max=5000), 
                precipitation=runif(40,min=2,max=20))

#Convert to data.table
require("data.table")
d <- data.table(d)

##Create your desired variables
r <- d[,list(altitude=mean(altitude,na.rm=T),
             max.ndvi=max(ndvi,na.rm=T),
             max.month=month[ndvi==max(ndvi,na.rm=T)],
             max.timestamp=timestamp[ndvi==max(ndvi,na.rm=T)],
             max.precipitation=precipitation[ndvi==max(ndvi,na.rm=T)],
             cv=sd(precipitation,na.rm=T)/mean(precipitation,na.rm=T)
),by=c("station","landcover","year")]

##Create lagged variables
r[, c(paste("cv.lag", 1:10, sep="")) := lapply(1:10, function(i) c(rep(NA, i), head(cv, -i))), by=list(station,landcover)]

#Create list of the cv variable positions plus station,landcover, and year positions
pos <- grep("cv", names(r))
pos <- c(1:3,pos)

#Melt lagged variables
require("reshape2")
r.long <- melt(r[,pos,with=F],id.vars=c("station","landcover","year"),variable.name="cv",value.name="cv.val")

#Merge back on max ndvi
r.long <- merge(r.long,r[,list(station,landcover,year,max.ndvi)],by=c("station","landcover","year"),all.x=T)

#Only use landcovers and lags with at least 3 non-missing observations
r.ref <- r.long[,list(Obs=sum(is.na(cv.val)==0)),by=c("landcover","cv")][Obs>2]
r.long <- merge(r.ref,r.long,by=c("landcover","cv"))

#Print out anovas
r.long[,print(summary(aov(max.ndvi~cv.val))),by=c("landcover","cv")]

#If you just care about pvalues, use this code
lmp <- function (modelobject) {
  if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
  f <- summary(modelobject)$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
}

#Print out results
r.long[,list(PVAL=lmp(lm(max.ndvi~cv.val))),by=c("landcover","cv")]