如何计算广义线性混合效应模型的中位绝对偏差(MAD)

Question

我知道我的问题与统计数据有关，但我正在R中寻找一个解决方案，所以我认为它适合这样做。

我建立了一个广义线性混合效应模型(GLMM)，该模型使用glmer函数从R中的lme4包建立了水产养殖站点周围物种丰富度的模型，该模型基于Zuur等人的重要解释变量。(2009年) https://rads.stackoverflow.com/amzn/click/com/1441927646。模式是：

Mod1 <- glmer(Richness ~ Distance + Depth + Substrate + Beggiatoa + 
        Distance*Beggiatoa + (1|Site/transect), family = poisson, data = mydata)

现在，我在不同的站点上收集了完整的数据集，我想评估这个模型在新数据集上的表现。

继简历上的问题之后，有人建议在新的数据集上寻找中位数绝对偏差(mad)。我尝试了来自stats包的R函数，但是我得到了以下错误消息：

Error in x[!is.na(x)] : object of type 'S4' is not subsettable
In addition: Warning messages:
1: In is.na(x) : is.na() applied to non-(list or vector) of type 'S4'
2: In is.na(x) : is.na() applied to non-(list or vector) of type 'S4'

，有人知道这里出了什么问题吗？难道mad stats 中的不能计算出GLMMs吗？如果是这样的话，是否有另一个包可以从GLMMs计算mad？

编辑：

为了让您了解我的数据，这里是dput(head(mydata))的输出，还请注意新的数据集中没有“基板”类别，"S“指的是”丰富性“：

structure(list(S = c(0, 1, 2, 3, 3, 2), Site = structure(c(1L, 
1L, 1L, 1L, 1L, 1L), .Label = c("BC", "BH", "GC", "IS", "Ref"
), class = "factor"), Transect = structure(c(4L, 4L, 4L, 4L, 
4L, 4L), .Label = c("10GC", "10IS", "10N", "10S", "11IS", "12IS", 
"13E", "1GC", "1N", "1W", "2E", "2GC", "2IS", "2N", "2W", "2WA", 
"3E", "3GC", "3IS", "3N", "3S", "4E", "4GC", "4IS", "4S", "4W", 
"5GC", "5IS", "5S", "6GC", "6IS", "6N", "6S", "6W", "7E", "7GC", 
"7IS", "8GC", "8IS", "8W", "9E", "9GC", "9IS", "9N", "RefBC1", 
"RefBC10", "RefBC11", "RefBC12", "RefBC2", "RefBC3", "RefBC4", 
"RefBC5", "RefBC6", "RefBC7", "RefBC8", "RefBC9", "X1", "X2"), class = "factor"), 
Distance = c(2, 20, 40, 80, 120, 160), Depth = c(40L, 40L, 
50L, 40L, 40L, 40L), Beggiatoa = c(2, 1, 1, 0, 0, 0)), .Names = c("S", 
"Site", "Transect", "Distance", "Depth", "Beggiatoa"), row.names = c(NA, 
6L), class = "data.frame")

Answer 1

对于样本内的误差，中位绝对偏差的计算将是

mad(residuals(fitted_model))

..。您可能需要residuals(fitted_model,type="response")，因为默认情况下，residuals会给出偏差残值(参见?residuals.merMod)。

如果您想查看样本外的错误，可以这样做：

pred <- predict(fitted_model,
                newdata = newdf,
                type = "response",
                re.form=~0)
mad(pred, center=newdf$S)

(re.form=~0指定要省略预测中的随机效果，这是您唯一的选择，除非您在也获得培训数据的站点/横断面上进行预测)