Очевидно, что это уже давно, но за последние несколько лет были достигнуты довольно большие успехи в работе с автокоррелированными данными.Я делаю то же самое в настоящее время, только биномиальный И с ошибкой missclassificaiton (идти большой или идти домой, верно?).Во всяком случае, я использовал код Скотта-Хавайра из пакета MRSea и Pirotta et al.2011 для определения значимых предикторов в автокоррелированных временных данных.Ниже я подхожу к этому.Я написал пользовательскую функцию подбора модели, которая медленно (чтобы избежать злоупотреблений и неправильного использования) идентифицирует значимые предикторы (с помощью обратного выбора QIC), а также функцию, которая выполняет повторные тесты Вальдса для определения значимости.В этом случае тип леса был действительно значительным, а виды - нет.
Для структуры блокировки, которую я использовал, вид, вероятно, более консервативен, чем вам нужно (типа леса может быть достаточно), но я склонен быть консервативным в своем выборе.Надеюсь, что это поможет, и, пожалуйста, всем цитируйте Pirotta et al 2011 и пакет MRSea (доступный на веб-сайте CREEM), если вы используете этот код.
Также обратите внимание, что функция anova может не подходить.Пиротта использовал функцию 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)