2012-02-07 4 views
3

나는 시계열이 여러 개 있는데, 각 시계열은 식물 종을 나타냅니다. 우디 밀도에 따라 패턴이 있다고 생각합니다. 높은 우디 밀도 종은 비가 오는 기간에 꽃이 피 웁니다. 낮은 기간 동안 우디 밀도가 높은 종입니다.여러 시계열의 차이를 테스트하는 방법 R

많은 종의 시계열과 우디 밀도 측정치를 사용하여이 패턴을 보여주기 위해 이것을 R으로 모델링하는 방법은 무엇입니까?

#Woody Density 
set.seed(69) 
wden<-round(c(rnorm(5,mean=50),rnorm(5,mean=90))) 
names(wden)<-c(paste("sp",1:10,sep="")) 
wden 

#Chuva 
rain<-c(150,100,50,40,20,20,30,50,70,100,150,200,150,100,50,30,20,20,40,50,70,100,150,200) 


#Flowering measures 
ydet<-c(10,10,10,10,20,40,50,40,20,10,10,10) 

#2 years for 5 low woody density and 5 high density species 
flowering<-matrix(NA,nrow=24, ncol=10,dimnames=list(paste("month",1:24,sep=""),paste("sp",1:10,sep=""))) 
for (i in 1:5) { 
       flowering[,i]<-round(c(ydet+rnorm(12,mean=5,sd=5),ydet+rnorm(12,mean=5,sd=5)),digits=2) 
       } 
for (i in 6:10) { 
       flowering[,i]<-round(c(rnorm(12,mean=30,sd=5),rnorm(12,mean=30,sd=5)),digits=2) 
       } 
#Changing objects to Time series 
flowering<-ts(flowering) 
#Plot series 
plot(flowering) 

#Making colors for wood density 
cores<-heat.colors(10,alpha=1) 
matplot(c(1:24),flowering,type="l",lwd=2,lty=1,xlab="Time",ylab="Flowering",col=cores[order(wden)]) 

#Plotting Rain Together with time series 
bargraph<-barplot(rain/max(rain)*100,xlab="Time",ylab="Rain") 
matlines(bargraph,flowering,type="l",lwd=2,lty=1,xlab="Time",ylab="Flowering",col=cores[order(wden)]) 
axis(1,at=bargraph,labels=1:24) 
axis(4,at=seq(0,100,by=10)) 
+0

'brota'와'dmad' 객체를 정의하는 것을 잊었을 것 같습니다. – nograpes

+0

작전, 코드 포르투갈어, 그때 더 많은 의미를 만들기 위해 일을 번역하려고했는데, 내가 모든 라인에서 이름을 변경하는 것을 잊지 것, 죄송합니다, 지금 바로 수정합니다 :) –

+0

글쎄, 나는 모든 것이 작동한다고 생각합니다. 감사 –

답변

2

실제로 당신이 http://stats.stackexchange.com에 시도하거나 [email protected] 메일 링리스트에 제안 될 수 있습니다 여기에

데이터가 같은 모습의 예입니다. 벌레가 조금 있습니다. 근본적인 문제는 두 시간 계열의 연관성이 인과 관계라는 것을 증명하는 것이 어렵다는 것입니다. (특히 둘 다 시간이 지남에 따라 규칙적으로 변동하는 경우) 동일한 기간에 단순히 변동하기 쉽고 따라서 줄을서는 것처럼 보일 수 있습니다 이것의 고전적인 예는 흑점주기와 뉴욕 증권 거래소와 같이 전혀 관련이없는 다양한 시계열입니다.

클래식 접근 방식은 각각의 시계열을 구별 할 수 없을 때까지 시계열 각각을 "희게"하는 것입니다 (주기적인 스플라인 또는 사인파 모델에 맞출 수도 있고 계절적 평균과 패턴의 차이를 따를 수도 있습니다) 백색 잡음을 측정 한 다음 시간차 (교차점 0, 즉 정규 상관 관계 또는 선도/후속 패턴을 나타내는 다른 차선)에서 상호 상관을 조사하십시오. 당신의 경우에는 상호 상관이 나무 밀도와 어떻게 다른지를 볼 수 있습니다.

"비오는 계절"과 "건기"로 데이터를 일괄 처리하고 더 표준적인 분석을 수행하면 (덩어리로 상관 관계를 대부분 제거 할 수 있음), (1) 데이터를보고이를 수행하는 것보다 선험적으로 계절을 나누십시오. (2) 이런 식으로 약간의 힘 및/또는 재미있는 미세 패턴을 잃을 수 있습니다. (3) 기본적인 유인 문제 중 일부는 여전히 존재한다 - 비와 관련된 꽃이 피고, 아니면 단지 비가 내리는 꽃이있다 계절?

좋은 예입니다.

0

분명히 이것은 기한이 지난 것입니다. 그러나 지난 몇 년 동안 자기 상관 데이터 작업에 상당한 진보가있었습니다. 나는 현재 같은 일을하고 있으며, 단지 2 항항과 missclassificaiton 오류 (큰 것 또는 집으로 돌아갈 것인가?)를 가지고있다. 어쨌든 MRSea Package와 Pirotta 외의 Scott-Hawyard 코드를 사용했습니다. 자기 상관 된 시간 데이터에서 중요한 예측자를 결정하기 위해 2011 년 아래는 제가 어떻게 접근하고 있는지입니다. 필자는 의미있는 예측 자 (역방향 QIC 선택을 통해)를 식별하기 위해 느린 (남용 및 오용을 방지하는 데 도움이되는) 사용자 정의 모델 피팅 함수를 작성했으며, 중요성을 판별하기 위해 Walds 테스트를 반복하는 피팅 함수를 작성했습니다. 이 경우 숲 유형은 실제로 중요하고 종은 그렇지 않습니다.

필자는 종을 사용해 왔지만, 필요한 것보다 더 보수적 일 수 있습니다 (포리스트 유형이면 충분할 수 있음). 그러나 나는 자신의 선택에 보수적 인 경향이 있습니다. 이 코드가 도움이 되었으면합니다.이 코드를 사용하면 Pirotta 외 2011 및 MRSea 패키지 (CREEM 웹 사이트에서 제공)를 인용하십시오.

또한 anova 기능이 적절하지 않을 수 있습니다. Pirotta는 geepack에서 anova.gee 함수를 사용했지만 그 함수는 새로운 버전에서 사라진 것으로 보이며 아직 대체품을 찾지 못했습니다. 제안을 환영합니다.

### Functions 

# This code calculates do backwards stepwise QIC to select the best model, 
# then repeated Walds tests to retain meaningful variables and CalcAUC to caluclate 
# the area under the ROC curve for each fitted model. All code based on 
# Pirotta E, Matthiopoulos J, MacKenzie M, Scott-Hayward L, Rendell L Modelling sperm whale habitat preference: a novel approach combining transect and follow data 
library(geepack) 
library(splines) 
library(AUC) 
library(PresenceAbsence) 
library(ROCR)   # to build the ROC curve 
library(mvtnorm)   # for rmvnorm used in predictions/plotting 


SelectModel=function(ModelFull){ 

    # Calculate the QIC of the full model 
    fullmodQ=QIC(ModelFull) 
    newQIC=0 
    terms=attr(ModelFull$terms,"term.labels") 


    while(newQIC != fullmodQ & length(terms)>1){ 


    # get all the terms for the full model 
    terms <- attr(ModelFull$terms,"term.labels") 
    n=length(terms) 

    newmodel=list() 
    newQIC=list() 

    newmodel[[1]]=ModelFull 
    newQIC[[1]]=fullmodQ 

    # Make n models with selection 
    for (ii in 1:n){ 
     dropvar=terms[ii] 
     newTerms <- terms[-match(dropvar,terms)] 
     newform <- as.formula(paste(".~.-",dropvar)) 
     newmodel[[ii+1]] <- update(ModelFull,newform) 
     newQIC[[ii+1]] =QIC(newmodel[[ii]]) 

    } 

    # Get the model with the lowest QIC 
    LowestMod=which.min(unlist(newQIC)) 

    if (LowestMod != 1){ 
     ModelFull=newmodel[[LowestMod]] 
     newQIC=min(unlist(newQIC)) 
    } else { 
     ModelFull=ModelFull 
     newQIC=min(unlist(newQIC)) 
    } 


    #end the model selection 


    } 
    return(ModelFull) 

} 

DropVarsWalds=function(ModelFull){ 

    # If no terms included return 
    if (length(attr(ModelFull$terms,"term.labels"))<2){ 
    NewModel='No Covariates to select from' 

    }else{ 


    OldModel=ModelFull 
    # Get the anova values 
    temp=anova(ModelFull) 

    # Make n models with selection 
    while(length(which(temp$`P(>|Chi|)`>.05))>0 & is.data.frame(temp)){ 


     # get the maximum value 
     dropvar=rownames(temp)[which.max(temp$`P(>|Chi|)`)] 

     # new formula for the full model 
     newform <- as.formula(paste(".~.-",dropvar)) 

     # new full model 
     ModelFull= update(ModelFull,newform) 

     # Get the model covariate names 
     terms <- attr(ModelFull$terms,"term.labels") 

     # # Get the anova values 
     # temp=anova(ModelFull) 

     temp=tryCatch({anova(ModelFull)}, error=function(e){e}) 


    } 

    NewModel=ModelFull 
    } 

    return(NewModel) 
} 


# functions to produce partial plots (largely based on MRSea, Scott-Hawyard) 
partialDF=function(mod, data, Variable){ 

    coefpos <- c(1, grep(Variable, colnames(model.matrix(mod)))) 
    xvals <- data[, which(names(data) == Variable)] 
    newX <- seq(min(xvals), max(xvals), length = 500) 
    eval(parse(text = paste(Variable, "<- newX", 
          sep = ""))) 
    response <- rep(1, 500) 
    newBasis <- eval(parse(text = labels(terms(mod))[grep(Variable, 
                 labels(terms(mod)))])) 
    partialfit <- cbind(rep(1, 500), newBasis) %*% coef(mod)[coefpos] 
    rcoefs <- NULL 
    try(rcoefs <- rmvnorm(1000, coef(mod), summary(mod)$cov.scaled), 
     silent = T) 
    if (is.null(rcoefs) || length(which(is.na(rcoefs) == 
             T)) > 0) { 
    rcoefs <- rmvnorm(1000, coef(mod), as.matrix(nearPD(summary(mod)$cov.scaled)$mat)) 
    } 
    rpreds <- cbind(rep(1, 500), newBasis) %*% t(rcoefs[, 
                 coefpos]) 
    quant.func <- function(x) { 
    quantile(x, probs = c(0.025, 0.975)) 
    } 
    cis <- t(apply(rpreds, 1, quant.func)) 

    partialfit <- mod$family$linkinv(partialfit) 
    cis <- mod$family$linkinv(cis) 

    fitdf=data.frame(x=newX, y=partialfit, LCI=cis[,1], UCI=cis[,2]) 
    colnames(fitdf)[1]=Variable 
    return(fitdf) 
} 

partialdf_factor=function(mod, data, variable){ 
    coeffac <- c(1,grep(variable, colnames(model.matrix(mod)))) 
    coefradial <- c(grep("LocalRadialFunction", colnames(model.matrix(mod)))) 
    coefpos <- coeffac[which(is.na(match(coeffac, coefradial)))] 
    xvals <- data[, which(names(data) == variable)] 
    newX <- sort(unique(xvals)) 
    newX <- newX[2:length(newX)] 
    partialfit <- coef(mod)[c(coefpos)] 
    rcoefs <- NULL 
    try(rcoefs <- rmvnorm(1000, coef(mod), summary(mod)$cov.scaled), 
     silent = T) 
    if (is.null(rcoefs) || length(which(is.na(rcoefs) == T)) > 0) { 
    rcoefs <- rmvnorm(1000, coef(mod), as.matrix(nearPD(summary(mod)$cov.scaled)$mat)) 
    } 

    if((length(coefpos))>1){ 

    rpreds <- as.data.frame(rcoefs[, c(coefpos)]) 
    fitdf_year=data.frame(vals=rpreds[,1]) 
    fitdf_year$FactorVariable=as.factor(paste(variable, levels(xvals)[1], sep = '')) 

    ############################ 
    # Recompile for plotting # 
    ############################# 


    for(jj in 2:ncol(rpreds)){ 
     temp=data.frame(vals=rpreds[,jj]) 
     temp$FactorVariable=as.factor(colnames(rpreds)[jj]) 
     fitdf_year=rbind(fitdf_year, temp) 
     rm(temp) 
    } 

    }else{ 

    rpreds <- rcoefs[,coefpos] 
    fitdf_year=data.frame(vals=rpreds) 
    fitdf_year$FactorVariable=colnames(model.matrix(mod))[coefpos[2:length(coefpos)]] 

    } 

    fitdf=fitdf_year 
    colnames(fitdf)[2]=variable 

    return(fitdf) 


} 


#Reshape your data 

temp=data.frame(flowering) 
flowering=data.frame(y=temp[,1], spp=colnames(temp)[1]) 

for(ii in 2:ncol(temp)){ 

    flowering=rbind(flowering, data.frame(y=temp[,ii], spp=colnames(temp)[ii])) 

} 

flowering$rain=rep(rain, ncol(temp)) 
flowering$woods=as.factor(c(rep('dense',nrow(temp)*5), rep('light',nrow(temp)*5))) 

fulmod=geeglm(formula=y~bs(rain)+woods, 
       family = gaussian, 
       id=spp, 
       data=flowering, 
       corstr = 'ar1') 

# Use stepwise QIC to select the best model 
QICmod=SelectModel(fulmod) 

# Use repeated walds to select signficinat variables 
sig_mod=DropVarsWalds(ModelFull = QICmod) 
관련 문제