setup
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')
实际答案
添加到我之前的答案，因为托马斯似乎对人们如何处理感兴趣factor以及如何使用引导程序获得置信区间。
处理因素
首先处理因子并不比处理数值变量困难多少。这里的区别在于
- 当绘制对数值变量的影响时，因子应设置为其基本水平（例如，对于
am作为一个因素，其值为 1)
- 绘制因子时，将所有数值变量设置为平均值，并将所有其他因子设置为基准水平。
获取因子基准水平的一种方法是factor[which.min(factor)]还有一个是factor(levels(factor)[0], levels(factor)). The ggeffects包使用一些与此类似的方法。
引导
现在，实践中的引导有容易的，也有困难的。可以使用参数、半参数或非参数引导程序。
非参数引导是最容易解释的。只需对原始数据集进行采样（例如 2/3、3/4 或 4/5。较少的数据可用于“良好”的较大数据集），使用该样本重新拟合模型，然后预测该新模型。然后重复该过程 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不是显示混合效果的数据集，因此引导置信区间will完全不正常（值的范围太分散）。另请注意，对于这样一个不稳定的数据集，相当多的引导程序无法收敛到合理的结果。
对于参数引导，我们可以转向lme4::bootMer。该函数需要一个merMod model (glmer or 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 i use geom_point and geom_errorbar我在哪里使用geom_line and 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)')
