r - 随着时间的推移从 gam 模型计算年利率？

Question

问题移至 CrossValidated

---

我试图表达 gam 中两个类别之间“增长速度”的差异造型。我的数据表示随着时间的推移 [0-100%] 的累积值，但我希望(为了与其他研究的可比性)以年度值来表示它们。我想知道它是否有意义，如果我的模型是 gam ，因此是非线性的？

以前，根据复合年增长率 (CARG) 从 gam 的起始值和终止值/年数计算每两组的“年度值”预测值 here .现在，我想知道是否有可能以某种方式从模型本身获取这些值(value)，并将其转化为年增长率？诸如“glm”中的“Betas”之类的东西表明一个预测变量比另一个更重要？

我的数据是这样的:我有两组，随着时间的推移森林损失的累积率(从没有森林损失到完全森林损失的范围在 0-100% 之间)。例如。集团yellow预测从 2000 年到 2017 年将失去约 40% 的森林，小组 blue已预测同时丢失 ~ 10%。我的虚拟数据可用 here在dput下面的表格。我的伪代码模型:forest_loss ~ year + location ，其中不同的组代表不同的位置。

# Restructure the data first
  dd$location = as.factor(dd$location)  # claim grouping variable as factor
  dd$forest_loss <- dd$forest_loss/100  # the the values between 0-1

# Plot the gam model
   ggplot(dd, aes(x = year, 
                y = forest_loss, 
                group = location,
               color = location)) +
  geom_smooth(method = 'gam', 
              method.args = list(family = "betar"),
              formula = y~s(x, bs = 'cs'))  + 
  geom_point()

我的虚拟模型的公式类似于:m1<-bam(y ~ s(x, by = grp) + grp, dd, family = betar) (我使用的是用于真实数据的表格，虽然我在那里有更多的预测变量)，基于 bam/gam .家庭betar是因为预测值的范围在 0-100(0-1)之间。

这是虚拟示例，因此部分代码 (gam) 可能会产生错误。但我的主要问题是如何从累积值本身来表达每组的年度增长？它是否已经可以从summary(m1)获得，有意义吗？

谢谢你的想法。

虚拟数据 (shared as suggested here using the dput )

   dd<- structure(list(year = c(2000L, 2001L, 2002L, 2003L, 2004L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 
2015L, 2016L, 2017L, 2000L, 2001L, 2002L, 2003L, 2004L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 
2015L, 2016L, 2017L), forest_loss = c(2, 2.5, 2.7, 2.9, 2.9, 3, 3.1, 3.1, 
3.1, 4, 5, 5.6, 5.8, 5.9, 6.7, 7.2, 9, 10, 3, 3.2, 3.4, 3.5, 
3.5, 7, 7.2, 8, 15, 19, 22, 25, 27, 29, 32, 33, 35, 40), location = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L), .Label = c("blue", "yellow"), class = "factor")), row.names = c(NA, 
-36L), class = "data.frame")

---

回复评论:

不知何故，tidymv::plot_smooths(m1)提供错误:

Error in `parse_expr()`:
! `x` must contain exactly 1 expression, not 0.
Run `rlang::last_error()` to see where the error occurred.

summary(m1) 的输出

Family: Beta regression(75.174) 
Link function: logit 

Formula:
y ~ s(x, by = grp)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -2.671      0.083  -32.19   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                 edf Ref.df Chi.sq  p-value    
s(x):grpblue   1.000  1.000  11.64 0.000647 ***
s(x):grpyellow 3.684  4.547 283.12  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.952   Deviance explained =   85%
-REML = -64.262  Scale est. = 1         n = 36

Gratia 的部分效果:gratia::draw(m1, scales = 'free')

Answer 1

如果您想使用 tidymv 中的 plot_smooth，您应该像这样指定 series 和 comparison 参数:

library(mgcv)
library(tidymv)
dd$location = as.factor(dd$location)  # claim grouping variable as factor
dd$forest_loss <- dd$forest_loss/100  # the the values between 0-1

m1<-bam(forest_loss ~ s(year, by = location) + location, dd, family = betar)
tidymv::plot_smooths(m1, year, location)

^{创建于 2022-09-09 reprex v2.0.2}