机器算法验证 - 在 R 中使用平方根链接模拟 GLM - 吾爱随笔录

在 R 中使用平方根链接模拟 GLM

机器算法验证 r 广义线性模型模拟链接功能

2022-04-14 12:42:47

我正在尝试使用基本函数来模拟拟合的 GLM，而不是使用被广泛质疑和回答的模拟（）和预测（）函数。当我将我的数学函数与模拟（）和预测（）函数进行比较时，我得到了不同的结果。很可能我做错了什么，但我似乎找不到错误。

首先，我生成带有倾斜因变量的随机相关数据：

library(MASS)
n <- 1500
d <- mvrnorm(n=n, mu=c(0,0,0,0), Sigma=matrix(.7, nrow=4, ncol=4) + diag(4)*.3)
d[,1] <- qgamma(p=pnorm(q=d[,1]), shape=2, rate=2)

接下来，我使用平方根链接拟合 GLM 以调整偏斜数据：

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4], family=gaussian(link="sqrt"))

接下来，我首先在线性尺度上生成预测值，然后在包含随机不确定性的逆尺度上生成预测值（最终我想使用源数据以外的数据作为输入值）：

p_lin <- m$coef[1] + m$coef[2]*d[,2] + m$coef[3]*d[,3] + m$coef[4]*d[,4]
p <- rnorm(n=n, mean=p_lin^2, sd=sd(p_lin^2 - d[,1]))

我将结果与模拟（）和预测（）函数进行比较：

par(mfrow=c(1,1), mar=c(4,2,2,1), pch=16, cex=0.8, pty="s")
xylim <- c(min(c(d[,1], p)), max(c(d[,1], p)))
plot(p, d[,1], xlab="predicted values", ylab="original data", xlim=xylim, ylim=xylim, col=rgb(0,0,0,alpha=0.1))
points(simulate(m)$sim_1, d[,1], col=rgb(0,1,0,alpha=0.1))
points(predict(m, type="response"), d[,1], col=rgb(1,0,0,alpha=0.1))
abline(a=0, b=1, col="red")

我的公式预测值出了什么问题？我在哪里可以找到预测（）和模拟（）函数中使用的数学和 R 表达式？是否有解释 GLM 模拟的链接，包括在 R 中应用的各种族/链接组合中的随机不确定性（最终我的下一步也是参数不确定性）。我找到了 GLM 模拟的一个很好的来源，但没有回答我的具体问题：http：/ /www.sagepub.com/upm-data/57233_Chapter_6.pdf

2个回答

我想出了一些问题的答案。

关于预测和模拟的数学表达式，可以通过以下代码获得（感谢 W. van der Elst 的提示）：

getS3method(c("predict"), class = "glm")
getS3method(c("simulate"), class = "lm")

关于predict功能，我错误地使用了 option type=”response”。这不包括我的目标的随机不确定性。它是在反变换线性尺度上的预测。这可以通过下图进行测试（在问题中说明的初始代码之后运行）：

plot(predict(m, type="response"), p_lin^2)

关于simulate函数，如果我从最初的问题看情节似乎是正确的：

plot(p, d[,1], xlab="predicted values", ylab="original data", xlim=xylim, ylim=xylim, col=rgb(0,0,0,alpha=0.1))
points(simulate(m)$sim_1, d[,1], col=rgb(0,1,0,alpha=0.1))

但是，如果我放大因变量：

d[,1] <- d[,1]*10000

并重新计算预测：

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4], family=gaussian(link="sqrt"))
p_lin <- m$coef[1] + m$coef[2]*d[,2] + m$coef[3]*d[,3] + m$coef[4]*d[,4]
p <- rnorm(n=n, mean=p_lin^2, sd=sd(p_lin^2 - d[,1]))

并绘制结果：

par(mfrow=c(1,1), mar=c(4,2,2,1), pch=16, cex=0.8, pty="s")
xylim <- c(min(c(d[,1], p)), max(c(d[,1], p)))
plot(p, d[,1], xlab="predicted values", ylab="original data", xlim=xylim, ylim=xylim, col=rgb(0,0,0,alpha=0.1))
points(simulate(m)$sim_1, d[,1], col=rgb(0,1,0,alpha=0.1))

我看到了不同的预测。

这似乎归因于我的通用公式方法和simulate()函数方法中的 sd 值之间的差异。回想一下通用公式方法中 sd 的计算：

sd=sd(p_lin^2 - d[,1])

在模拟函数 R 中计算 sd 的方式与一般公式方法（对于高斯族）不同：

vars <- deviance(m)/df.residual(m)
if (!is.null(m$weights)) vars <- vars/m$weights # the m$weights seems similar to the m$fitted.values multiplied by about 4
fitted(m) + rnorm(n, sd = sqrt(vars))

我不明白为什么要计算 sd 然后除以m$weights. 为什么 sd 是向量而不是单个值？simulate()帮助文本指出： “由 lm 或 glm(family = "gaussian") 拟合的线性模型的方法假设已提供的任何权重都与误差方差成反比。” 我似乎无法理解这句话的意思。如果我使用该simulate()函数运行多个模拟，模拟看起来非常相似：

plot(simulate(m)$sim_1,simulate(m,2)$sim_2)

如果我不使用 sqrt-link 函数，我得到的模拟似乎更好地反映了随机不确定性，因为当我运行它们两次时它们不太相同：

library(MASS)
rm(list=ls())
set.seed(2)
n <- 1500
d <- mvrnorm(n=n, mu=c(0,0,0,0), Sigma=matrix(.7, nrow=4, ncol=4) + diag(4)*.3)
d[,1] <- qgamma(p=pnorm(q=d[,1]), shape=2, rate=2)
d[,1] <- d[,1]*10000
m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4])
plot(simulate(m)$sim_1,simulate(m,2)$sim_2)

什么方法是正确的？有什么区别？（我应该将此作为单独的问题发布吗？）

有什么区别？

该型号

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4], family=gaussian(link="sqrt"))

是

y_{i} \sim N (μ_{i}, σ^{2}), \sqrt{μ_{i}} = {\vec{β}}^{⊤} {\vec{x}}_{i}

$y_i \sim N(\mu_i,\sigma^2), \quad \sqrt{\mu_i} = \vec{\beta}^\top\vec{x}_i$

它是用迭代加权最小二乘法估计的。因此你会发现使用了更多的迭代

library(MASS)
set.seed(2)
n <- 1500
d <- mvrnorm(n=n, mu=c(0,0,0,0), Sigma=matrix(.7, nrow=4, ncol=4) + diag(4)*.3)
d[,1] <- qgamma(p=pnorm(q=d[,1]), shape=2, rate=2)

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4], family=gaussian(link="sqrt"))
m$iter
#R> [1] 5

它是一个广义线性模型，weights元素包含收敛时的工作权重。另一方面这个模型

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4])

是

y_{i} \sim N (μ_{i}, σ^{2}), μ_{i} = {\vec{β}}^{⊤} {\vec{x}}_{i}

$y_i \sim N(\mu_i,\sigma^2), \quad \mu_i = \vec{\beta}^\top\vec{x}_i$

并在一次迭代中用最小二乘法估计。

m1 <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4])
m1$iter # takes an extra iteration in `glm.fit` loop
#R> [1] 2

m2 <- lm(formula=d[,1] ~ d[,2] + d[,3] + d[,4])
stopifnot(all.equal(coef(m1), coef(m2))) # gives the same

您开始使用的模型也不是吗？

关于

我不明白为什么要计算 sd 然后除以m$weights. 为什么 sd 是向量而不是单个值？

那么以下

线性模型的方法由lm或glm(family = "gaussian")假设已提供的任何权重与误差方差成反比。

是有道理的，例如，如果你有一个平均结果 $n_i$ 具有相同协变量的观测值。平均响应的方差将与 $n_i$ 您将作为重量提供。但是，当您使用其他链接函数和恒等函数时，它不仅提供权重，还提供收敛的工作权重。我不确定这是否有意义。

其它你可能感兴趣的问题

上一篇用于变化检测的功率鞅：M 变为零？下一篇总和差异的置信区间