显然，这早就该这么做了，但在过去几年中，在处理自相关数据方面取得了一些相当大的进步。我目前也在做同样的事情，只有二项式 AND 错误分类错误（变大还是回家，对吗？）。无论如何，我使用了来自 MRSea 包和 Pirotta 等人的 Scott-Hawyard 的代码。 2011 年确定自相关时间数据中的重要预测因子。以下是我的处理方法。我编写了一个自定义模型拟合函数，该函数很慢（有助于避免滥用和误用）来识别有意义的预测变量（通过向后 QIC 选择）以及进行重复的 Walds 测试以确定显着性的函数。在这种情况下，森林类型确实很重要，而物种则不重要。
对于我使用物种的阻塞结构，它可能比你需要的更保守（森林类型可能就足够了），但我在选择上倾向于保守。希望这有帮助，如果您使用此代码，请大家引用 Pirotta et al 2011 和 MRSea 包（可在 CREEM 网站上获得）。
另请注意，方差分析函数可能不合适。 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)