如何将 gam() 合并到 Lattice 包中的 xyplot() 中？答案

【问题标题】：How to incorporate gam() in xyplot() from Lattice package?如何将 gam() 合并到 Lattice 包中的 xyplot() 中？
【发布时间】：2014-06-01 08:15:55
【问题描述】：

我正在尝试将来自mgcv 包的广义加法模型gam() 合并到R 中lattice 包中的xyplot() 函数或coplot() 函数。

可以通过选择臭氧数据在http://statweb.stanford.edu/~tibs/ElemStatLearn/ 中找到数据。

这是我的内核平滑代码。

    ozonedata=ozone 
    Temperature=equal.count(ozonedata$temperature,4,1/2)
    Wind=equal.count(ozonedata$wind,4,1/2)

    xyplot(ozone^(1/3) ~ radiation | Temperature * Wind, data = ozonedata, as.table = TRUE, 
           panel = function(x, y, ...) {panel.xyplot(x, y, ...);panel.loess(x, y)}, 
           pch = 20,xlab = "Solar Radiation", ylab = "Ozone (ppb)")

或

    coplot((ozone^(1/3))~radiation|temperature*wind,data=ozonedata,number=c(4,4),
           panel = function(x, y, ...) panel.smooth(x, y, span = .8, ...),
           xlab="Solar radiation (langleys)", ylab="Ozone (cube root ppb)")

广义加性模型生成如下。

    gam_ozone = gam(ozone^(1/3)~s(radiation)+s(temperature)+s(wind),data=ozonedata,method="GCV.Cp")

现在我无法将 gam() 的拟合组合到格子图中。

【问题讨论】：

标签： r lattice gam mgcv

【解决方案1】：

我认为这应该可行。请注意，我们将gam_ozone 对象作为参数传递给xyplot，称为gam。这将允许我们在面板功能中访问它。

xyplot(ozone^(1/3) ~ radiation | Temperature * Wind, data = ozonedata, 
    as.table = TRUE, gam=gam_ozone, 
    panel = function(x, y, gam,...) {

        xx <- dimnames(trellis.last.object())
        ww <- arrayInd(packet.number(), .dim=sapply(xx, length))
        avgtmp <- mean(xx[[1]][[ww[1]]])
        avgwind <- mean(xx[[2]][[ww[2]]])

        gx<-seq(min(x, na.rm=T), max(x, na.rm=T), length.out=50)
        gy<-predict(gam, data.frame(radiation=gx, 
            temperature=avgtmp, wind=avgwind))

        panel.xyplot(x, y, ...);
        panel.xyplot(gx, gy, ..., col="red", type="l");
    }, 
    pch = 20,xlab = "Solar Radiation", ylab = "Ozone (ppb)"
)

现在为了使用gam 进行预测，您必须为每个面板找到一个用于风和温度的值。我决定做的只是为每个瓦片范围取中间值。所以我使用Deepayan-approved, undocumented feature 来获取每个当前带状疱疹的范围，并调用dimnames。然后我找到带有packet.number() 的当前面板一旦我有了范围，我就会采取手段获得平均值。

我不会使用我们传入的gam 模型来预测每个面板的曲线。我根据观察值计算了 50 个x 值的范围，然后从gam 预测一个新行。

最后我用panel.xyplot 绘制原始数据，然后绘制gam 预测线。

【讨论】：