如何将多元 GAM 模型的 y 轴从平滑值更改为实际值？

Question

我正在使用多元 GAM 模型来详细了解多个地区的雾趋势。雾是由低于某个阈值（

但是，我现在面临的挑战是，我真的希望 y 轴成为实际的能见度观测值，而不是居中的平滑值。 有趣的是，看到协变量相对于该位置的平均能见度如何影响能见度，但是对于平均能见度不同的多个位置（因此 0 点在提高或降低可见度的意义不大）。

为了比较多个位置的结果，我正在尝试对 y 轴进行实际能见度观察，然后在我们有兴趣查看的能见度阈值（400 m）处画一条线 评估低于该阈值的预测变量值是什么样的（例如，什么温度与低于 400 m 的能见度相关）。

总的来说，在 GAM 和 R 方面，我仍然是初学者，但到目前为止我已经找到了一些有用的部分。

到目前为止有用的东西：

尝试1.如何为模型中的每个变量提取gam fit Extracting data used to make a smooth plot in mgcv

尝试2.如何使用predict函数重构单变量模型 http://zevross.com/blog/2014/09/15/recreate-the-gam-partial-regression-smooth-plots-from-r-package-mgcv-with-a-little-style/

尝试 3. 如何使用“拟合”获得一些看起来像可见性观察的 y 轴 - 虽然我不认为这是 正确的方法，因为我没有考虑拦截 http://gsp.humboldt.edu/OLM/R/05_03_GAM.html

模拟数据

install.packages("mgcv") #for gam package
require(mgcv)
install.packages("pspline")
require(pspline)


#simulated GAM data for example
dataSet <- gamSim(eg=1,n=400,dist="normal",scale=2)
visibility <- dataSet[[1]]
temperature <- dataSet[[2]]
dewpoint <- dataSet[[3]]
windspeed <- dataSet[[4]]


#Univariable GAM model
gamobj <- gam(visibility ~  s(dewpoint))
plot(gamobj, scale=0, page=1, shade = TRUE, all.terms=TRUE, cex.axis=1.5, cex.lab=1.5, main="Univariable Model: Dew Point")
summary(gamobj)
AIC(gamobj)
abline(h=0)

露点单变量模型 https://imgur.com/1uzP34F

ATTEMPT 2 -- 使用单变量模型预测函数，但不改变 y 轴

#dummy var that spans length of original covariate
maxDP <-max(dewpoint)
minDP <-min(dewpoint)
DPtrial.seq <-seq(minDP,maxDP,length=3071)
DPtrial.seq <-data.frame(dewpoint=DPtrial.seq)

#predict only the DP term 
preds <- predict(gamobj, type="terms", newdata=DPtrial.seq, se.fit=TRUE)

#determine confidence intervals
DPplot <-DPtrial.seq$dewpoint
fit <-preds$fit
fit.up95 <-fit-1.96*preds$se.fit
fit.low95 <-fit+1.96*preds$se.fit

#plot
plot(DPplot, fit, lwd=3,
 main="Reconstructed Dew Point Covariate Plot")

#plot confident intervals
polygon(c(DPplot, rev(DPplot)), 
    c(fit.low95,rev(fit.up95)), col="grey",
    border=NA)

lines(DPplot, fit,  lwd=2)
rug(dewpoint) 

重构的露点协变量图 https://imgur.com/VS8QEcp

ATTEMPT 3 -- 使用“拟合”更改 y 轴，但不考虑截距

plot(dewpoint,fitted(gamobj), main="Fitted Response of Y (Visibility) Plotted Against Dew Point")
abline(h=mean(visibility))
rug(dewpoint)

Y 对露点绘制的拟合响应 https://imgur.com/RO0q6Vw

最终，我想要一条水平线，我可以在其中研究相对于 400 米的预测变量，而不仅仅是响应变量的平均值。这样，它将在平均可见度不同的多个站点之间具有可比性。最重要的是，它需要用于多个协变量！

Gavin Simpson 在几篇文章中解释了该方法，但不幸的是，我真的不明白在使用预测函数时如何保持其他协变量的平均值不变：

Changing the Y axis of default plot.gam graphs

对此方法的任何更深入的解释都会非常有帮助！！！

Answer 1

我不确定这会有多大帮助，因为您的 Q 比我们通常在 SO 上想要的更开放，但是，就这样吧。

首先，我认为考虑对响应变量进行建模会有所帮助，我认为这是当前的可见性。这将是一个连续变量，以 0 为界（也许数据永远不会达到零？），这表明将数据建模为条件分布

• gamma (family = Gamma(link = 'log')) 的可见性永远不会为零。
• Tweedie (family = tw()) 用于包含零的数据。

另一种方法是对雾的发生进行建模；如果这被定义为能见度 family = binomial() 将数据建模为条件分布的伯努利。

确定了建模方法后，我们需要对响应进行建模。这应该使用多重回归类型的方法来完成，其中 GAM 包括多个预测变量。通过这种方式，您可以估计每个潜在预测变量对响应的影响，同时控制其他预测变量的影响。如果您一次只使用一个预测变量（例如dewpoint）执行此操作，那么该变量可以很好地“解释”可能由于另一个预测变量 windspeed 导致的数据变化，而您不会知道它。

此外，预测变量之间很可能存在交互，如果它们存在，您需要对其进行控制，这只能在

中完成

然后，为了最终解决问题的症结，在拟合多预测模型以“解释”可见性后，您需要从模型中预测可能的条件集。要查看在其他预测变量有影响的模型中可见性如何随dewpoint 变化，您需要将其他变量固定在一些合理的值；一种选择是将它们设置为它们的平均值（或在任何因子预测变量的情况下为模态值），或指示该变量的典型值的其他值。为此，您必须使用您的领域知识。

如果您在模型中有交互，那么您需要改变交互中的两个变量，同时将所有其他变量固定在某些值。

假设您没有交互并且对dewpoint 感兴趣，但该模型还包括windspeed。用于拟合模型的值的平均风速可以从拟合模型的cmX 组件中找到。您可以根据观察到的windpseed 值计算此值，或将其设置为您想要使用的某个已知数字。用m 表示拟合，用df 表示其中包含您的数据的数据框，然后我们可以创建新数据以在dewpoint 的范围内进行预测，同时保持windspeed 不变。

mn.windspd <- m$cmX['windspeed']
## or
mn.windspd <- with(df, mean(windspeed))
## or set it some some value
mn.windspd <- 10 # say

那你就可以了

preddata <- with(df,
                 expand.grid(dewpoint = seq(min(dewpoint),
                                            max(dewpoint),
                                            length = 300),
                             windspeed = mn.windspd))

然后你使用它从拟合模型中进行预测：

pred <- predict(m, newdata = preddata, type = "link", se.fit = TRUE)
pred <- as.data.frame(pred)

现在我们想要将这些预测重新放到响应尺度上，并且我们想要一个置信区间，因此我们必须在反向转换之前先创建它：

ilink <- family(m)$linkinv
pred <- transform(pred,
                  Fitted = ilink(fit),
                  Upper  = ilink(fit + (2 * se.fit)),
                  Lower  = ilink(fit - (2 * se.fit)),
                  dewpoint = preddata = dewpoint)

现在您可以将dewpoint 对响应的影响可视化，同时保持windspeed 不变。

在您的情况下，您必须将其扩展为保持 temperature 不变，但以相同的方式完成

mn.windspd <- m$cmX['windspeed']
mn.temp <- m$cmX['temperature']
preddata <- with(df,
                 expand.grid(dewpoint = seq(min(dewpoint),
                                            max(dewpoint),
                                            length = 300),
                             windspeed = mn.windspd,
                             temperature = mn.temp))

然后按照上面的步骤进行预测。

对于一个或两个不同的变量，我的 gratia 包中有一个函数 data_slice()，它将为您完成上述 expand.grid() 的工作，因此您不必指定平均值其他协变量：

preddata <- data_slice(m, 'dewpoint', n = 300)

从技术上讲，这会在数据中找到最接近中值的值（对于不变的协变量）。如果你想要手段，那就去做

fixdf <- data.frame(windspeed = mn.windspd, temperature = mn.temp)
preddata <- data_slice(m, 'dewpoint', data = fixdf, n = 300)

如果你有交互，比如dewpoint 和windspeed，那么你需要改变两个变量。再次使用expand.grid()，这很容易：

mn.temp <- m$cmX['temperature']
preddata <- with(df,
                 expand.grid(dewpoint = seq(min(dewpoint),
                                            max(dewpoint),
                                            length = 100),
                             windspeed = seq(min(windspeed),
                                             max(windspeed),
                                             length = 300),
                             temperature = mn.temp))

这将创建一个 100 x 100 的协变量值网格来预测，同时保持温度恒定。

对于data_slice()，您需要这样做：

fixdf <- data.frame(temperature = mn.temp)
preddata <- data_slice(m, 'dewpoint', 'windpseed',
                       data = fixdf, n = 300)

并将其扩展到您想要改变的更多协变量，使用expand.grid() 遵循此模式也很容易；我还没有实现超过 2 个在 data_slice 中变化的变量。