将具有截断功率基的 GAM 表示为混合模型

机器算法验证混合模式广义加法模型样条

2022-03-24 00:45:34

截断的幂基础如下所示：

y = θ_{0} + θ_{1} x + θ_{2} x^{2} + . . . + θ_{d} x^{d} + \sum_{k = 1}^{K} θ_{d k} (x - κ_{k})_{+}^{d}

$y = \theta_0 + \theta_1x + \theta_2x^2 + ... + \theta_dx^d + \sum_{k=1}^K \theta_{dk}(x-\kappa_k)_+^d$

显然，您可以很容易地将系数分成固定和随机效应，如下所示：

y = \overset{fixed effects}{\overset{⏞}{θ_{0} + θ_{1} x + θ_{2} x^{2} + . . . + θ_{d} x^{d}}} + \overset{random effects}{\overset{⏞}{\sum_{k = 1}^{K} θ_{d k} (x - κ_{k})_{+}^{d}}}

$y = \overbrace{\theta_0 + \theta_1x + \theta_2x^2 + ... + \theta_dx^d}^{\text{fixed effects}} + \overbrace{\sum_{k=1}^K \theta_{dk}(x-\kappa_k)_+^d}^{\text{random effects}}$

我在这方面找到了很多资料，但都认为这种分离是给定的，没有人解释为什么会这样。为什么必须以这种方式分离固定效应和随机效应，在这种情况下随机效应分布是什么？

1个回答

您在问题中讨论的模型可以写成

y = X β + F b + e

$y = X \beta + F b+e$ 在哪里

X

$X$ 是列等于的矩阵

1, x, x^{2}, x^{3}, . . .

$1, x, x^2, x^3,...$ 和

F

$F$ 是一个矩阵，其列是通过计算截断的多项式获得的。

那么（惩罚）目标函数是：

Q_{p} = ‖ y - X β + F b ‖^{2} + k ‖ b ‖^{2}

$Q_{p} = \|y - X \beta + F b\|^2 + k\|b\|^{2}$ 只有 s 系数被缩小。

b

$b$

为了计算和，我们需要求解以下惩罚正态方程系统：你可以比较 eqs 的系统。上面有一个，例如，这里https://en.wikipedia.org/wiki/Mixed_model（估计会话）。方差分量是和和。 $\beta$ $b$

[\begin{array}{lll} X^{'} X & X^{'} F \\ F^{'} X & F^{'} F + k I \end{array}] [\begin{array}{ll} β \\ b \end{array}] = [\begin{array}{ll} X^{'} y \\ F^{'} y \end{array}]

$\left[ \begin{array}{lll} X'X & X'F \\ F' X & F'F + kI \end{array} \right] \left[ \begin{array}{ll} \beta\\ b \end{array} \right] = \left[ \begin{array}{ll} X'y \\ F'y \end{array} \right]$

σ^{2} = v a r (e)

$\sigma^2 = var(e)$

τ^{2} = v a r (b)

$\tau^2 = var(b)$

k = σ^{2} / τ^{2}

$k = \sigma^{2}/\tau^{2}$

为什么必须以这种方式分离固定效应和随机效应：您会注意到，在 Henderson 的混合模型方程中，随机效应也被“惩罚”（项）。 $G^{-1}$

这种情况下的随机效应分布是什么：我们假设和 $b \sim N(0, \tau^{2} I)$ $e \sim N(0, \sigma^{2} I)$

我希望我的回答能有所帮助，并且我的符号是正确的。

编辑

评论：为什么tpf部分需要被惩罚？

像往常一样，惩罚控制平滑度和数据拟合之间的权衡（见下图，其中我用 15 个 2 度 TPF 基和不同级别的 k 参数对相同的数据进行平滑处理）。这适用于所有惩罚平滑技术。

我们为什么要做这一切？

使混合效应模型表示法方便的是这样一个事实，即模型（包括最佳平滑量）可以使用标准 lmm 例程计算（下面我使用 nlme ...请注意，我假设您有一个计算 tpf_bases 的函数）。

# Simulate some data
n  = 30
x  = seq(-0, 2*pi, len = n)
ys = 2 * sin(x)
y  = rnorm(n, ys, 0.5)

# Create bases
Bs = tpf_bases(x, ndx = 10, deg = 2) 
X  = Bs$X
Z  = Bs$Z

# Organize for lme
dat     = data.frame(X1 = X[, 2], X2 = X[, 3], y = y)
dat$Z   = Z
dat$all = (1:n) * 0 + 1

# Fit lme
fit = lme(y ~ X1 + X2, random = list(all = pdIdent( ~ Z - 1)), data = 
dat)

# Extract coefficients & get fit
beta.hat  = fit$coef$fixed
b.hat     = unlist(fit$coef$random)
f.hat     = X %*% beta.hat + Z %*% b.hat

# Plot results
plot(x, y, main = "LME-based optimal fit")
lines(x, f.hat, col = 'red')

其它你可能感兴趣的问题

上一篇t 分布是指数族的成员吗？下一篇为什么低概率事件很少发生的非重言式解释？