问如何使用R检验几个时间序列之间的差异

EN

### 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)