我正在gam从mgcv包中使用模型并存储结果,model到目前为止,我一直在使用平滑组件plot(model).我最近开始使用ggplot2并喜欢它的输出.所以我想知道,是否可以使用ggplot2绘制这些图表？
这是一个例子:
x1 = rnorm(1000)
x2 = rnorm(1000)
n = rpois(1000, exp(x1) + x2^2)
model = gam(n ~ s(x1, k=10) + s(x2, k=20), family="poisson")
plot(model, rug=FALSE, select=1)
plot(model, rug=FALSE, select=2)
我感兴趣s(x1, k=10)而s(x2, k=20)不是合适.
部分答案:
我更深入地挖掘plot.gam和mgcv:::plot.mgcv.smooth和建立了自己的功能,提取从平滑部件的预测效果和标准误差.它没有处理所有选项和案例,plot.gam因此我只将其视为部分解决方案,但它对我来说效果很好.
EvaluateSmooths = function(model, select=NULL, x=NULL, n=100) {
if (is.null(select)) {
select = 1:length(model$smooth)
}
do.call(rbind, lapply(select, function(i) {
smooth = model$smooth[[i]]
data = model$model
if (is.null(x)) { …在使用mgcv包运行gam模型时,我遇到了一条奇怪的错误消息,我无法理解:
"model.frame.default中的错误(公式=死亡~pm10 +滞后(resid1,1)+:变量长度不同(找到'Lag(resid1,1)')".
模型1中使用的观察数量与偏差残差的长度完全相同,因此我认为此误差与数据大小或长度的差异无关.
我在网上找到了一个相当有关的错误信息在这里,但后没有得到充分的答案,因此它是不利于我的问题.
可重复的示例和数据如下:
library(quantmod)
library(mgcv)
require(dlnm)
df <- chicagoNMMAPS
df1 <- df[,c("date","dow","death","temp","pm10")]
df1$trend<-seq(dim(df1)[1]) ### Create a time trend
model1<-gam(death ~ pm10 + s(trend,k=14*7)+ s(temp,k=5),
data=df1, na.action=na.omit, family=poisson)
resid1 <- residuals(model1,type="deviance")
model1_1 <- update(model1,.~.+ Lag(resid1,1), na.action=na.omit)
model1_2<-gam(death ~ pm10 + s(trend,k=14*7)+ s(temp,k=5) + Lag(resid1,1), data=df1,
na.action=na.omit, family=poisson)
这两个模型都产生了相同的错误消息.
我正在使用GAM来模拟逻辑回归中的时间趋势.然而,我想从中提取拟合样条曲线,将其添加到另一个模型中,不能用于GAM或GAMM.
因此我有两个问题:
随着时间的推移,我怎样才能更顺畅,以便在让模特找到其他结点的同时强迫一个结处于特定位置？
如何从拟合的GAM中提取矩阵,以便我可以将其用作不同模型的推算？
我正在运行的模型类型如下:
gam <- gam(mortality.under.2~ maternal_age_c+ I(maternal_age_c^2)+
s(birth_year,by=wealth2) + wealth2 + sex +
residence + maternal_educ + birth_order,
data=colombia2, family="binomial")
我已经阅读了GAM的大量文档,但我还不确定.任何建议都非常感谢.
我正在使用 mgcv 包编写 GAM，该包使用实地考察期间获得的数据和从 Sentinel 卫星拍摄的图像来预测岛上两种不同物种的洞穴丰度和分布。调查了101个地块。在 66 个地块中记录了属于物种 1 的 922 个洞穴，在 8 个地块中记录了属于物种 2 的洞穴29 个。
我对物种 1 使用了负二项式分布，因为使用泊松分布会导致模型过度分散。最大模型是：
gam(Species_1 ~ s(x, y, bs="ts") +
Sentinel2_band_1 + Sentinel2_band_2 + Sentinel2_band_3 + Sentinel2_band_4 + Sentinel2_band_5 +
Sentinel2_band_6 + Sentinel2_band_7 + Sentinel2_band_8 + Sentinel2_band_9 + Sentinel2_band_10 +
I(Sentinel2_band_1^2) + I(Sentinel2_band_2^2) + I(Sentinel2_band_3^2) + I(Sentinel2_band_4^2) + I(Sentinel2_band_5^2) +
I(Sentinel2_band_6^2) + I(Sentinel2_band_7^2) + I(Sentinel2_band_8^2) + I(Sentinel2_band_9^2) + I(Sentinel2_band_10^2) +
aspect + elevation + slope +
I(aspect^2) + I(elevation^2) + I(slope^2) +
aspect:elevation + …我尝试将游戏适合我拥有的一些数据框。全部减去一项工作。它失败并出现错误：
Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : A term has fewer unique covariate combinations than specified maximum degrees of freedom
我在互联网上查了一下，但无法真正弄清楚到底出了什么问题。我所有的 7 个以上数据帧都运行没有问题。然后我跑了 epiR::epi.cp(srtm[-c(1,7,8)])，它给了我这个输出：
$cov.pattern
id n curv_plan curv_prof dem slope ca
1 1 1 1.113192e-02 3.991046e-03 3909 43.601479 5.225853
2 2 1 -2.686749e-03 3.474989e-03 3312 35.022511 4.418310
3 3 1 -1.033450e-02 -4.626922e-03 3326 36.678623 4.421465
4 4 1 -5.439283e-03 2.066148e-03 4069 31.501045 3.887526
5 5 1 -2.602015e-03 -1.249511e-04 3021 37.199219 5.010560
6 6 1 1.068216e-03 1.216902e-03 2844 44.694374 …我想在两个约束，以适应GAM模型数据simultatenously：（1）拟合单调（增加），（2）配合经过一个固定的点，比如说，(x0,y0)。
到目前为止，我设法让这两个约束分开工作：
对于 (1)，基于mgcv::pcls() 文档示例，通过使用mgcv::mono.con()来获得足以满足单调性的线性约束，并通过mgcv::pcls()使用约束来估计模型系数。
对于（2），基于这篇文章，通过使用模型公式中的偏移项将节点位置 x0 处的样条值设置为 0 +。
但是，我很难同时结合这两个约束。我想一种方法是mgcv::pcls()，但我既不能解决 (a) 使用偏移将节点位置 x0 处的样条值设置为 0 + 的类似技巧，也不能 (b) 设置相等约束（我认为）可以产生我的（2）约束设置）。
我还注意到，对于我的约束条件 (2)，将结点位置 x0 处的样条值设置为 0 的方法产生了奇怪的摆动结果（与不受约束的 GAM 拟合相比）——如下所示。
模拟一些数据
library(mgcv)
set.seed(1)
x <- sort(runif(100) * 4 - 1)
f <- exp(4*x)/(1+exp(4*x))
y <- f + rnorm(100) * 0.1
dat <- data.frame(x=x, y=y)
GAM 无约束（用于比较）
k <- 13
fit0 …我有一个gam我知道可以正常使用的模型R,但是当我尝试train使用该caret包" "相同的模型时,它会返回一个错误,指出输入数据列是列表.有谁理解这个？
我正在运行的代码如下:
library("caret")
library("mgcv")
a <- gam(RW ~ s(Temp0.grd) + s(mld.grd) + s(mean_depth.grd) +
s(land_dist.grd) + s(slope.grd) + s(npp.grd),
data=mydata,
family=binomial)
all.data.gam.train <-
train(form=RW ~ s(Temp0.grd) + s(mld.grd) + s(mean_depth.grd) +
s(land_dist.grd) + s(slope.grd) + s(npp.grd),
data=mydata,
method='gam',
family=binomial
)
第一个gam模型工作正常,但是train返回以下错误:
Error in model.frame.default(form = RW ~ s(Temp0.grd) + s(mld.grd) + s(mean_depth.grd) + :
invalid type (list) for variable 's(Temp0.grd)'
直接在公式上运行model.frame.default也会产生这个错误,因此严格来说问题不在于火车.
mydata看起来如下:
> class(mydata)
[1] "data.frame"
> class(mydata$Temp0.grd)
[1] "numeric"
> …使用mgcv的惩罚样条,我希望在示例数据中获得10 /年的有效自由度(EDF)(整个期间为60).
library(mgcv)
library(dlnm)
df <- chicagoNMMAPS
df1<-subset(df, as.Date(date) >= '1995-01-01')
mod1 <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow)
,family=quasipoisson,na.action=na.omit,data=df1)
在示例数据中,由edf测量的时间的基本维度为56.117,小于每年10.
summary(mod1)
Approximate significance of smooth terms:
edf Ref.df F p-value
s(time) 56.117 67.187 5.369 <2e-16 ***
s(temp) 2.564 3.204 0.998 0.393
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.277 Deviance explained = 28.2%
GCV score = 1.1297 Scale est. = 1.0959 n = 2192
手动我将通过提供如下的平滑参数来更改edf a
mod1$sp
s(time) …我想gam在mgcv包中使用函数:
x <- seq(0,60, len =600)
y <- seq(0,1, len=600)
prova <- gam(y ~ s(x, bs='cr')
我可以设置结数s()吗？然后我可以知道花键使用的结点在哪里？谢谢!
我对根据全国抽样调查数据进行的GAM回归很有趣。我对此文章感兴趣 。我选择了感兴趣的变量生成DF：
nhanesAnalysis <- nhanesDemo %>%
select(fpl,
age,
gender,
persWeight,
psu,
strata)
然后，据我所知，我使用以下代码生成了加权DF：
library(survey)
nhanesDesign <- svydesign( id = ~psu,
strata = ~strata,
weights = ~persWeight,
nest = TRUE,
data = nhanesAnalysis)
假设我只选择具有age?30：
ageDesign <- subset(nhanesDesign, age >= 30)
现在，我将使用拟合GAM模型（fpl ~ s(age) + gender）mgcv package。是否可以通过weights参数或使用svydesignobject来实现ageDesign？
编辑
我想知道从svyglm对象推断计算的权重并将其用作weightsGAM中的参数是否正确。