R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据
  MFqKcJmbshKy 2023年11月19日 48 0

我们使用广义线性模型(Generalized Linear Models,简称GLM)来研究客户的非正态数据,并探索非线性关系。GLM是一种灵活的统计模型,适用于各种数据类型和分布,包括二项分布、泊松分布和负二项分布等非正态分布。通过GLM,我们可以对非正态数据进行建模和预测,并且能够处理计数数据,如客户购买数量、网站点击次数等。GLM还允许引入自变量的非线性效应,从而更好地拟合与响应变量之间的复杂关系。这使得GLM成为处理非正态数据和非线性关系的强大工具。

泊松回归和伽马回归 - 探索联系

如果我们查看火车与机动车碰撞数据,我们会发现一个有趣的模式。

 

 

library(readr)
......

train <- read_csv("s_agresti.csv")

train_plot <- ggplot(train,
   ......

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布

仅仅通过观察,我们就可以看出方差随预测变量而变化。此外,我们处理的是计数数据,它具有自己的分布,即泊松分布。然而,如果我们坚持使用lm进行分析会怎样呢?

 

 

train_lm <-......odel(train_lm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_02

预测值和观测值之间不匹配。部分原因是这里的响应变量在残差中不是正态分布的,而是泊松分布,因为它是计数数据。

泊松回归

具有泊松误差的广义线性模型通常具有对数链接,尽管也可以具有恒等链接。例如,

 

 

pois_tib <- tibble(x = rep(0:40,2),
            ......
  geom_col(position = position_dodge())

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_03

上面显示了两个泊松分布,一个均值为5,另一个均值为20。请注意它们的方差如何变化。

对数链接(例如ŷ=ea+bx̂=eβ+αx)是一个自然的拟合方法,因为它不能得到小于0的值。因此,在这种情况下,我们可以这样做:

 

 

train_glm <- glm(collisions ~ km_trtravel, 
           ......

然后,我们可以重新评估模型的假设,包括过分离。请注意,下面的QQ图并没有什么实际意义,因为这不是正态分布。

 

 

check_model(train_glm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_04

那么...残差怎么办呢?鉴于残差不是正态分布的,使用qqnorm图几乎没有意义。拟合残差关系仍然可能看起来很奇怪。

使用广义线性模型的分位数残差

评估广义线性模型(以及许多其他模型形式)的一种方法是查看其分位数残差。 因此,首先让我们使用DHARMa生成一些模拟残差。

 

 

res <- simulateResiduals(train_glm)

我们可以绘制这些图表,并进行非参数拟合检验。

 

 

plotQQunif(res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_05

很好,拟合效果不错。忽略异常值测试,因为在更详细的观察中我们发现没有异常值。

我们还可以查看预测与量化残差图。

 

 

plotResiduals(res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_06

 

 

check_overdispersion(train_glm) |> plot()

 

 

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_07

查看glm结果

让我们来看看模型结果。

 

 

summary(train_glm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_08

注意,在这里我们看到了标准的glm输出,我们可以像处理任何对数变换一样解释系数。我们还有一个离散参数,描述了均值和方差之间的关系。对于泊松分布,它的值为1。

最后,我们可以绘图。

 

 

train_plot +
......"log")))

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_09

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_10

Gamma回归

然而,我们的数据通常不是离散的。考虑来自Zuur的蛤蜊数据。

 

 

clams <- read_delim("ams.txt", delim = "\t") %>%
  mutate(MONTH = factor(MONTH))

AFD(无灰干质量)与月份和长度有什么关系?显然这里存在非线性关系。

 

 

clam_plot <- ggplot(clams,
      ......

clam_plot

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_11

现在,看起来我们应该用对数变换的模型进行拟合,但是…

 

 

clam_lm <- lm(log(......

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_12

显然存在明显的问题。即使对AFD取对数后的qq图也不好,残差拟合图也不好。Gamma glm采用其逆函数作为其规范连接,但它们通常也可以使用对数连接。

 

 

clam_gamma <- glm(AFD ~ ......
                  data = clams)

 

 

check_model(clam_gamma)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_13

还有

 

 

clam_res <- simulateR......res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_14

 

 

ploals(clam_res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_15

好的,也许不是很好。但这主要是由于高值的稀疏性导致的,所以没关系。

我们可以使用predict进行绘图,在这里分别绘制每个月的图。

 

 

clam_plot +......
  facet_wrap(~MONTH)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_16

我们还可以查看其他属性。

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_17

 

 

summary(clam_gamma)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_18

我们可以重新参数化伽马分布,使得均值=形状/速率。在这种情况下,我们使用该均值和形状参数化伽马分布。离散参数是1/形状。

但是,为了更容易理解,伽马的方差随均值的平方成比例地扩展。离散参数越大,方差扩展得越快。

最后,我们可以使用纳吉尔克计的伪R2来计算R2。

 

 

# fit
r2(clam_gamma)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_19

这是正态的吗?

你可能会问为什么这里使用伽马分布而不是正态分布?我们可以用正态误差和对数链接进行glm拟合。

 

 

clam_glm_norm <- glm(AFD ......
                     data = clams)

一种判断的方法是寻找过离散。

 

 

norm_res <- simulateRe......orm_res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_20

 

 

plotuals(norm_res)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_21

我们可以看到QQ图很好。而且predobs也不糟糕(特别是与上面相比)。这是一些很好的证据,表明这里可能只需要正态误差和对数链接。

逻辑回归

让我们来看看我们的小鼠感染隐孢子虫的例子。请注意,数据被限制在0和1之间。

 

 

mouse <- read_csv...... Porportion)) +
  geom_point()

mouse_plot

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_22

这是因为虽然N是每个样本的总小鼠数量,但是我们不能有超过N的感染!实际上,每只老鼠就像一次抛硬币。它是否被感染了。

二项分布

二项分布有两个参数,成功的概率和硬币投掷的次数。得到的分布始终介于0和1之间。考虑使用不同概率进行15次硬币投掷的情况。

 

R

bin_tibble <- tibble(outcome = rep(0:15, 2),......
  geom_col(position = position_dodge())

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_23

我们也可以将x轴的范围调整为0到1,来表示比例。

或者,考虑相同的概率,但是不同次数的硬币投掷。

 

R

bin_tibble <- tibble(outcome = rep(0:15, 2),......
  geom_col(position = position_dodge())

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_24

你可以看到两个参数都会影响分布的形状。

二项式逻辑回归

在二项逻辑回归中,我们主要是估计获得正面的概率。然后我们以权重的形式提供(而不是估计)试验次数。这里使用的典型链接函数是logit函数,因为它描述了一个在0和1之间饱和的逻辑函数。

在R中,我们可以使用两种形式来参数化二项逻辑回归 - 这两种形式是等价的,因为它们将结果扩展为成功次数和总试验次数。

 

R

mouse_glm_cbind <- glm(cbind(Y,......
                 data = mouse)

第二种方式使用权重来表示试验次数。

 

R

mouse_glm <- glm(Porport......
                 data = mouse)

这两个模型是相同的。

从这一点开始,工作流程与以往一样 - 假设检验、分析和可视化。

 

R

checl(mouse_glm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_25

 

R

binduals(mouse_glm, ......

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_26

 

R

res_bin <- sim......

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_27

 

R

plotRes_bin)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_28

 

R

summary(moglm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_29

 

R

r2(mouse_glm)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_30

注意,离散参数为1,就像泊松分布一样。

 

R

ggplot(mouse,
     ......
              method.args = list(family = binomial))

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_31

Beta回归

最后,我们经常会遇到受限数据,但这些数据不是从二项式分布中抽取的 - 也就是说,并不存在独立的“硬币翻转”,如果你愿意的话。

考虑以下关于服用不同补充剂时锻炼后钠摄入比例的分析,2300是推荐摄入量,所以我们将其标准化为这个值。

 

R

sodium <- read_csv("laake.csv")

 

R

ggplot(sodium,
 ......
  geom_boxplot()

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_泊松分布_32

现在,让我们使用Beta回归来观察这个结果。

 

R

sodium_beta <- beta......
                       data = sodium)

 

R

soditmb <- glmmTMB(Porport......
                       data = sodium)
chec......a_tmb)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_33

 

R

plotQQunif(sodium_beta_tmb)

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_34

然后我们可以继续进行所有我们通常的测试和可视化。例如 -

 

R

emmeans(sodium_b......
  confint(adjust = "none")

R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_数据_35

如果我们有一个连续的协变量,我们可以获得拟合值和误差,并将它们放入模型中。


R语言非线性回归和广义线性模型:泊松回归、伽马回归、逻辑回归、Beta回归分析机动车事故、小鼠感染、蛤蜊数据、补剂锻炼钠摄入数据_拟合_36


【版权声明】本文内容来自摩杜云社区用户原创、第三方投稿、转载,内容版权归原作者所有。本网站的目的在于传递更多信息,不拥有版权,亦不承担相应法律责任。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@moduyun.com

  1. 分享:
最后一次编辑于 2023年11月19日 0

暂无评论

推荐阅读
MFqKcJmbshKy