Как проверить разницу между несколькими временными рядами, используя 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))
2 ответа
Я мог бы предложить вам попробовать это на http://stats.stackexchange.com/ или на r-sig-ecology@r-project.org список рассылки. Это немного банка червей. Основная проблема заключается в том, что трудно доказать, что связь двух временных рядов носит причинно-следственный характер, поскольку (особенно, когда оба регулярно меняются во времени) им легко просто колебаться примерно в один и тот же период и, следовательно, казаться выстраивающимися в линию (классические примеры этого - цикл солнечных пятен и различные совершенно не связанные временные ряды, такие как Нью-Йоркская фондовая биржа).
Классический подход к этому состоял бы в том, чтобы "отбелить" каждый из временных рядов (вы могли бы подогнать периодический сплайн или синусоидальную модель или просто взять отличия шаблона от сезонного среднего) независимо, пока каждый из них не будет неотличим от белого шума, затем изучите взаимные корреляции между временными рядами (при нулевом запаздывании, то есть регулярных корреляциях или при других запаздываниях, чтобы представить ведущую / следующую модель). В вашем случае вы, вероятно, захотите увидеть, как взаимные корреляции меняются в зависимости от плотности древесины.
В качестве альтернативы, вы можете просто объединить данные в "сезон дождей" и "сухой сезон" и сделать более стандартный анализ (вы избавитесь от большей части корреляции путем объединения), но (1) было бы неплохо иметь априори разделение времен года вместо того, чтобы делать это, просматривая данные; (2) таким образом вы можете потерять некоторую мощность и / или интересные мелкомасштабные шаблоны; (3) некоторые из основных проблем со здоровьем все еще существуют - цветение связано с дождем или просто с сезоном дождей?
Хороший пример, хотя.
Очевидно, это уже давно, но за последние несколько лет в работе с автокоррелированными данными были достигнуты довольно значительные успехи. Я делаю то же самое в настоящее время, только биномиальное И с ошибкой missclassificaiton (идти большой или идти домой, верно?). Во всяком случае, я использовал код Скотта-Хавярда из пакета MRSea и Pirotta et al. 2011 для определения значимых предикторов в автокоррелированных временных данных. Ниже я подхожу к этому. Я написал пользовательскую функцию подбора модели, которая медленна (чтобы избежать злоупотреблений и неправильного использования) для определения значимых предикторов (посредством обратного выбора QIC), а также функцию, которая выполняет повторные тесты Вальдса для определения значимости. В этом случае тип леса был действительно значительным, а виды - нет.
Для блокирующей структуры, которую я использовал, она, вероятно, более консервативна, чем вам нужно (тип леса может подойти), но я склонен быть консервативным в своем выборе. Надеюсь, что это поможет, и, пожалуйста, всем цитируйте Pirotta et al 2011 и пакет MRSea (доступный на веб-сайте CREEM), если вы используете этот код.
Также обратите внимание, что функция анова может не подходить. Пиротта использовал функцию anova.gee из пакета geepack, но эта функция, похоже, исчезла из более новых версий, и я пока не нашел замены. Предложения приветствуются.
### 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)