设置
library(lme4)
library(MuMIn)
database <- mtcars
database$wt <- scale(mtcars$wt)
database$am <- factor(mtcars$am) ## <== note the difference here. It is a factor not numeric
model.1 <- glmer(vs ~ wt + am + (1|carb), database, family = binomial, na.action = "na.fail")
model.1.set <- dredge(model.1, rank = "AICc")
top.models.1 <- get.models(model.1.set,subset = delta<10)
model.1.avg <- model.avg(top.models.1)
nPoints <- 100
wt_pred_data <- data.frame(wt = seq(min(database$wt), max(database$wt), length = nPoints),
am = database$am[which.min(database$am)], #Base level for the factor
var = 'wt')
am_pred_data <- data.frame(wt = mean(database$wt),
am = unique(database$am),
var = 'am')
pred_data <- rbind(wt_pred_data, am_pred_data)
rm(wt_pred_data, am_pred_data)
pred_data$vs <- predict(model.1.avg, newdata = pred_data, re.form = NA, type = 'response')
实际答案
添加到我之前的答案中，因为托马斯似乎对如何处理factors 以及如何使用引导程序获得置信区间感兴趣。
处理因素
首先处理因子并不比处理数值变量难多少。这里的区别在于
- 在绘制对数值变量的影响时，应将因子设置为它们的基本水平（例如，
am作为一个因子，这将是一个值 1）
- 绘制因子时，将所有数值变量设置为其平均值，将所有其他因子设置为其基本水平。
获得因子基本水平的一种方法是factor[which.min(factor)]，而另一种方法是factor(levels(factor)[0], levels(factor))。该ggeffects包使用一些类似于此的方法。
自举
现在，实践中的引导从容易到困难。可以使用参数、半参数或非参数引导程序。
非参数引导是最容易解释的。一个人只需简单地从原始数据集（比如 2/3、3/4 或 4/5。Less 可用于“好”的较大数据集）中提取样本，使用该样本重新拟合模型，然后预测该新模型。然后重复该过程 N 次，并将其用于估计标准偏差或分位数，并将其用于置信区间。似乎没有实现的方法MuMIn来为我们处理这个问题，所以我们似乎必须自己处理这个问题。
通常代码会变得非常混乱，使用函数可以使其更清晰。令我沮丧的是MuMIn然而，这似乎有问题，所以下面是一种非功能性的方法。在这段代码中，我选择了 4/5 的样本大小，因为数据集的大小相当小。
### ###
## Non-parametric bootstrapping ##
## Note: Gibberish with ##
## singular fit! ##
### ###
# 1) Create sub-sample from the dataset (eg 2/3, 3/4 or 4/5 of the original dataset)
# 2) refit the model using the new dataset and estimate model average using this dataset
# 3) estimate the predicted values using the refitted model
# 4) refit the model N times
nBoot <- 100
frac <- 4/5 #number of points in each sample. Better datasets can use less.
bootStraps <- vector('list', nBoot)
shutup <- function(x) #Useful helper function for making a function shut up
suppressMessages(suppressWarnings(force(x)))
ii <- seq_len(n <- nrow(database))
nn <- ceiling(frac * n)
nb <- nn * nBoot
samples <- sample(ii, nb, TRUE)
samples <- split(samples, (nn + seq_len(nb) - 1) %/% nn) #See unique((nn + seq_len(nb) - 1) %/% nn) # <= Gives 1 - 100.
#Not run:
# lengths(samples) # <== all of them are 26 long! ceiling(frac * n) = 26!
# Run the bootstraps
for(i in seq_len(nBoot)){
preds <- try({
# 1) Sample
d <- database[samples[[i]], ]
# 2) fit the model using the sample
bootFit <- shutup(glmer(vs ~ wt + am + (1|carb), d, family = binomial, na.action = "na.fail"))
bootAvg <- shutup(model.avg(get.models(dredge(bootFit, rank = 'AICc'), subset = delta < 10)))
# 3) predict the data using the new model
shutup(predict(bootAvg, newdata = pred_data, re.form = NA, type = 'response'))
}, silent = TRUE)
#save the predictions for later
if(!inherits(preds, 'try-error'))
bootStraps[[i]] <- preds
# repeat N times
}
# Number of failed bootStraps:
sum(failed <- sapply(bootStraps, is.null)) #For me 44, but will be different for different datasets, samples and seeds.
bootStraps <- bootStraps[which(!failed)]
alpha <- 0.05
# 4) use the predictions for calculating bootstrapped intervals
quantiles <- apply(do.call(rbind, bootStraps), 2, quantile, probs = c(alpha / 2, 0.5, 1 - alpha / 2))
pred_data[, c('lower', 'median', 'upper')] <- t(quantiles)
pred_data[, 'type'] <- 'non-parametric'
请注意，这当然完全是胡言乱语。拟合是奇异的，因为mtcars不是显示混合效应的数据集，因此自举置信区间将完全不合常理（值的范围过于分散）。还要注意，对于这样一个不稳定的数据集，相当多的引导程序无法收敛到合理的东西。
对于参数引导，我们可以转向lme4::bootMer. 此函数采用单个merMod模型（glmer或lmer结果）以及要在每个参数改装上评估的函数。所以创建这个函数bootMer可以处理剩下的事情。我们对预测值感兴趣，所以函数应该返回这些。注意功能的相似之处，与上述方法
### ###
## Parametric bootstraps ##
## Note: Singular fit ##
## makes this ##
## useless! ##
### ###
bootFun <- function(model){
preds <- try({
bootAvg <- shutup(model.avg(get.models(dredge(model, rank = 'AICc'), subset = delta < 10)))
shutup(predict(bootAvg, newdata = pred_data, re.form = NA, type = 'response'))
}, silent = FALSE)
if(!inherits(preds, 'try-error'))
return(preds)
return(rep(NA_real_, nrow(pred_data)))
}
boots <- bootMer(model.1, FUN = bootFun, nsim = 100, re.form = NA, type = 'parametric')
quantiles <- apply(boots$t, 2, quantile, probs = c(alpha / 2, 0.5, 1 - alpha / 2), na.rm = TRUE)
# Create data to be predicted with parametric bootstraps
pred_data_p <- pred_data
pred_data_p[, c('lower', 'median', 'upper')] <- t(quantiles)
pred_data_p[, 'type'] <- 'parametric'
pred_data <- rbind(pred_data, pred_data_p)
rm(pred_data_p)
再次注意，由于奇点，结果将是胡言乱语。在这种情况下，结果将过于确定，因为奇异性意味着模型在已知数据上过于准确。所以在实践中，这将使每个间隔的范围为 0 或足够接近以至于没有区别。
最后我们只需要绘制结果。我们可以facet_wrap用来比较参数和非参数的结果。再次注意，对于这个特定的数据集，比较完全无用的置信区间是非常无用的。
请注意，对于我使用的因子am和geom_point我geom_errorbar使用的位置geom_line以及geom_ribbon数值，与数值变量的连续性质相比，为了更好地表示因子的分组性质
#Finaly we can plot our result:
# wt
library(ggplot2)
ggplot(pred_data[pred_data$var == 'wt', ], aes(y = vs, x = wt)) +
geom_line() +
geom_ribbon(aes(ymax = upper, ymin = lower), alpha = 0.2) +
facet_wrap(. ~ type) +
ggtitle('gibberish numeric plot (caused by singularity in fit)')
# am
ggplot(pred_data[pred_data$var == 'am', ], aes(y = vs, x = am)) +
geom_point() +
geom_errorbar(aes(ymax = upper, ymin = lower)) +
facet_wrap(. ~ type) +
ggtitle('gibberish factor plot (caused by singularity in fit)')
