R에 대한 루프의 시계열 지연 실행 R

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 

나는 스테이션 당 연간 최대 NDVI에 강수량의 변동 계수의 지연 0시 10분의 효과를 발견하고 싶습니다. 는 기본적으로 내 코드는 다음과 같습니다

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 지연에는 모든 누락 값 (landcover level detail)이 있다고 생각합니다. 이 코드는 cv.lags (landcover 클래스 내에서)가 최소한 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")]