我正在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)
这两个模型都产生了相同的错误消息.
我正在尝试使用他们网站上的代码将Bioconductor安装到R中.当我输入代码时(见下文),我收到一条错误消息,说某些软件包无法更新,安装路径是不可写的.
> ## try http:// if https:// URLs are not supported
> source("https://bioconductor.org/biocLite.R")
Bioconductor version 3.4 (BiocInstaller 1.24.0), ?biocLite for help
> biocLite()
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.4 (BiocInstaller 1.24.0), R 3.3.2 (2016-10-31).
installation path not writeable, unable to update packages: Matrix, mgcv,
生存
我可以通过转包/安装包来安装这些包.
> utils:::menuInstallPkgs()
trying URL 'https://www.stats.bris.ac.uk/R/bin/windows/contrib/3.3/Matrix_1.2-8.zip'
Content type 'application/zip' length 2775038 bytes (2.6 MB)
downloaded 2.6 MB
trying URL 'https://www.stats.bris.ac.uk/R/bin/windows/contrib/3.3/mgcv_1.8- 16.zip'
Content type 'application/zip' length 2346257 bytes (2.2 MB)
downloaded 2.2 MB
trying URL 'https://www.stats.bris.ac.uk/R/bin/windows/contrib/3.3/survival_2.40-1.zip'
Content …几年前的这个主题描述了如何提取用于绘制拟合gam模型的平滑组件的数据.它有效,但只有当有一个平滑变量时才有效.我有多个平滑变量,不幸的是我只能从系列的最后一个中提取平滑.这是一个例子:
library(mgcv)
a = rnorm(100)
b = runif(100)
y = a*b/(a+b)
mod = gam(y~s(a)+s(b))
summary(mod)
plotData <- list()
trace(mgcv:::plot.gam, at=list(c(25,3,3,3)),
#this gets you to the location where plot.gam calls plot.mgcv.smooth (see ?trace)
#plot.mgcv.smooth is the function that does the actual plotting and
#we simply assign its main argument into the global workspace
#so we can work with it later.....
quote({
#browser()
plotData <<- c(plotData, pd[[i]])
}))
plot(mod,pages=1)
plotData
我试图让两个估计的平滑函数a和b,但列表plotData只给我估计b.我已经研究了plot.gam …
我正在使用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 …我对阅读gamm4包的源代码时注意到的一个事实感到困惑。
\n它在两个地方从mgcv导入内部函数。
\n一处位于函数中（此处gamm4.setup链接到代码）
G <- mgcv:::gam.setup(formula,pterms,\n data=data,knots=knots,sp=NULL,\n min.sp=NULL,H=NULL,absorb.cons=TRUE,sparse.cons=0,gamm.call=TRUE)\n另一个地方是 gamm4 函数（此处链接到代码）：
\n var.summary <- mgcv:::variable.summary(gp$pf,dl,nrow(mf)) ## summarize the input data\n在描述文件中，它在“取决于”下列出了 mgcv（在此处链接到完整的描述文件）：
\nDepends: R (>= 2.9.0), methods, Matrix, lme4 (>= 1.0), mgcv (>= 1.7-23)\n此外，在 NAMESPACE 中它导入了 mgcv，尽管我认为这与我的问题无关。此处链接到命名空间。
\n如果我在正在开发的包中执行完全相同的操作，则使用另一个包中的内部函数会导致警告，并留下以下警告R CMD check：
\xe2\x9d\xaf checking dependencies in R code ... WARNING\n Unexported objects imported by \':::\' calls:\n \xe2\x80\x98mgcv:::gam.setup\xe2\x80\x99 \xe2\x80\x98mgcv:::variable.summary\xe2\x80\x99\n …我在 ggplot2 中使用 stat_smooth 函数，决定我想要“拟合优度”，并为此使用了 mgvc gam。我突然想到我应该检查以确保它们是相同的模型（stat_smooth vs mgvc gam），所以我使用下面的代码进行检查。从表面上看，它们有不同的结果，正如情节所证明的那样（情节：stat_smoother gam (red), mgcv gam (black)）。但是，我不知道为什么他们有不同的结果。两者之间的某些默认参数是否不同？是不是 gam 是在数字 x 上运行而 stat_smooth 是在 POSIXct x 上运行的（如果是这样 - 我不知道该怎么办）？看起来 stat_smooth 更平滑，但 k 值是相同的...
我认为有几篇关于如何在 ggplot2 中绘制 gam 输出的帖子，但我真的很想知道为什么 stat_smooth 和 mgcv 首先给出不同的结果。我对 GAM（和 R）很陌生，所以很可能我错过了一些简单的东西。但是，我在问之前确实谷歌并搜索了这个论坛。
我的数据有点大，无法轻松共享，所以我使用了一个示例数据集 - 我已经将源代码放在代码中，以及dput()下面的所有内容，然后是我的sessionInfo()。
我试图提出一个质量问题，但这只是我的第二个问题。曾经。所以，建设性的批评是值得赞赏的。
谢谢！
library(readxl)
library(data.table)
library(ggplot2)
library(scales)
library(mgcv)
stackOF_data <- read_excel("mean-daily-flow-cumecs-vatnsdals.xlsx", sheet = "Data")
stackOF_data <- data.table(stackOF_data)
stackOF_data <- stackOF_data[,.(timeseries=as.POSIXct(Date,format("%Y-%m-%d")),mdf)]
a <- stackOF_data[,.(x=as.numeric(timeseries),y=mdf)]
a1 <- gam(y~s(x, k=100, bs="cs"),data=a)
a2=data.table(gam_mdf= …