机器算法验证 - 如果我希望模型截距具有全局平均值的解释，我应该使用哪种分类变量编码？ - 吾爱随笔录

如果我希望模型截距具有全局平均值的解释，我应该使用哪种分类变量编码？

机器算法验证混合模式固定效应模型分类编码

2022-03-23 12:17:04

在混合效应模型中

y_{i j} = β_{00} + β_{01} x_{1 i} + β_{02} x_{2 i} + β_{03} x_{3 i} + u_{i} + ϵ_{i j}

$y_{ij} = \beta_{00} + \beta_{01} x_{1i} + \beta_{02} x_{2i} + \beta_{03} x_{3i} + u_i + \epsilon_{ij}$

其中是虚拟变量，编码离散（多项式）变量的水平，它有两个以上的水平（这里是四个），我想给截距的解释（全局）总体平均值，即。 $x_1, x_2, x_3$ $\tilde{x}$ $\beta_{00}$ $E(y_{ij})$

现在协变量的编码方式是将截距解释为的参考类别的平均值。 $\tilde{x}$

有没有办法做到这一点？

我找到了关于效果编码的一个很好的概述，但这种类型的平均编码不是其中的一部分。

编辑：我只记得如何对只有两个类别的变量\ tilde {然后我们有模型 $\tilde{x}$

y_{i j} = β_{00} + β_{01} x_{1 i} + u_{i} + ϵ_{i j}

$y_{ij} = \beta_{00} + \beta_{01} x_{1i} + u_i + \epsilon_{ij}$

其中虚拟定义为如果并且它是如果，其中是与。 $x_{1i}$ $(1-p)$ $\tilde{x} =1$ $(-p)$ $\tilde{x} =0$ $p$ $\tilde{x}=1$

编辑 2：根据Robert Long 的回复，的每个级别的观察次数相同时，可以使用偏差编码。但是，我正在寻找一种可能具有不等类概率实现偏差编码的代码，证明这种编码不估计全局平均值。我怀疑需要对偏差编码假人进行某种类别加权（就像我对上面的两类案例所做的那样）。 $\tilde{x}$ $\tilde{x}$ $\tilde{x}$

# Code to assess deviation coding for multinomial $xt$    
library(MASS)
library(dplyr)
n = 1000
set.seed(13)
xt = rmultinom(n, 1, c(1/3, 1/3, 1/3))
xt = as.factor( apply( t(t(xt) * c(1,2,3)), 2, sum) )
X <- model.matrix(~ xt)
betas <- c(3, 1, 2)
Y <- X %*% betas + rnorm(n)
mean(Y)

lm(Y ~ xt) %>% coef()   # default treatment coding

contrasts(xt) <- contr.sum(3) # specify deviation coding
lm(Y ~ xt) %>% coef()

编辑 3：最初的问题是“如果我希望模型截距能够解释全局平均值，使用哪种效果编码（分类编码）？” 标题错误地暗示我的目标可以通过单独的效果编码来实现。Robert Long 的答案适用于平衡类别，然后应该使用偏差编码。

1个回答

如果数据是平衡的，那么偏差编码应该起作用。

让我们看一个简单的例子：

set.seed(13)
dt <- expand.grid(X1 = LETTERS[1:3], reps = 1:5)
X <- model.matrix(~ X1, dt)
betas <- c(3, 1, 2)
dt$Y <- X %*% betas + rnorm(nrow(dt))
mean(dt$Y)

[1] 4.11413

所以我们希望截距为 4.11413

如果我们用默认编码拟合模型，我们会得到：

lm(Y ~ X1, dt) %>% coef()   # default treatment coding

(Intercept)         X1B         X1C 
  3.3430627   0.2867999   2.0264018

但是现在如果我们使用偏差编码，我们会得到

contrasts(dt$X1) <- contr.sum(3) # specify deviation coding
lm(Y ~ X1, dt) %>% coef()

(Intercept)         X11         X12 
  4.1141299  -0.7710672  -0.4842673

如果数据不平衡，那么您将需要进行一些事后调整。

编辑：解决数据不平衡时该怎么做。

在这种情况下，使用默认处理编码比使用偏差编码更容易：

> set.seed(1)
> dt1 <- expand.grid(X1 = LETTERS[1:1], reps = 1:5)
> dt2 <- expand.grid(X1 = LETTERS[2:2], reps = 1:3)
> dt3 <- expand.grid(X1 = LETTERS[3:3], reps = 1:2)
> dt <- rbind(dt1, dt2, dt3)
> table(dt$X1)

A B C 
5 3 2

所以群体是不平衡的。

> X <- model.matrix(~ X1, dt)
> betas <- c(2, 3, 1)
> dt$Y <- 4 + X %*% betas + rnorm(nrow(dt), 0, 1)
> mean(dt$Y)

[1] 7.232203

所以我们想通过事后计算来恢复 7.23，这可以很容易地实现

> coef(lm(Y ~ X1, dt))[1] + betas[2] * table(dt$X1)[2]/nrow(dt) + betas[3] * table(dt$X1)[3]/nrow(dt)

(Intercept) 
   7.22927

请注意，由于组中的不平衡和随机误差的组合，结果并不准确。随着误差接近零，结果变得准确。即使有错误，结果也是无偏的，正如我们从蒙特卡罗模拟中看到的那样：

n.sim <- 1000
vec.sim <- numeric(n.sim)

for (i in 1:n.sim) {
  
  set.seed(i)

  dt$Y <- 4 + X %*% betas + rnorm(nrow(dt), 0, 1)

  vec.sim[i] <- mean(dt$Y) - (coef(lm(Y ~ X1, dt))[1] + betas[2] * table(dt$X1)[2]/nrow(dt) + betas[3] * table(dt$X1)[3]/nrow(dt))

}

hist(vec.sim)
mean(vec.sim)

[1] -0.003418483

编辑：如评论中所述，我们真的应该使用模型中的系数估计，这样做将使计算准确：

> coef(lm(Y ~ X1, dt))[1] + coef(lm(Y ~ X1, dt))[2] * table(dt$X1)[2]/nrow(dt) + coef(lm(Y ~ X1, dt))[3] * table(dt$X1)[3]/nrow(dt)
(Intercept) 
   7.232203

其它你可能感兴趣的问题

上一篇建议用于方差卡方检验的样本量下一篇分析部分交叉的设计