首页 文章

如何在lm中指定参数估计值之间的关系？

提问于
浏览 次
4

使用lm，我想拟合模型：y = b0 b1 * x1 b2 * x2 b1 * b2 * x1 * x2

我的问题是：如何指定相互作用的系数应该等于系数乘以主效应？

我已经看到将系数设置为特定值，您可以使用offset（）和I（），但我不知道如何指定系数之间的关系 .

这是一个简单的模拟数据集：

n <- 50 # Sample size
x1 <- rnorm(n, 1:n, 0.5) # Independent variable 1
x2 <- rnorm(n, 1:n, 0.5) # Independent variable 2
b0 <- 1 
b1 <- 0.5
b2 <- 0.2
y <- b0 + b1*x1 + b2*x2 + b1*b2*x1*x2 + rnorm(n,0,0.1)

为了拟合模型1：y = b0 b1 * x1 b2 * x2 b3 * x1 * x2，我会使用：

summary(lm(y~ x1 + x2 + x1:x2))

但是我如何拟合模型2：y = b0 b1 * x1 b2 * x2 b1 * b2 * x1 * x2？

两个模型之间的主要差异之一是要估计的参数数量 . 在模型1中，我们估计4个参数：b0（截距），b1（变量1的斜率），b2（变量2的斜率）和b3（变量1和2之间的相互作用的斜率） . 在模型2中，我们估计3个参数：b0（截距），b1（变量1的斜率和变量1和2之间相互作用的斜率的一部分）和b2（变量2的斜率和变化的斜率的一部分） . 变量之间的互动.1和2）

我想这样做的原因是，当调查x1和x2之间是否存在显着的相互作用时，模型2，y = b0 b1 * x1 b2 * x2 b1 * b2 * x1 * x2，可以是比y = b0 b1 * x1 b2 * x2 .

非常感谢！

玛丽

r lm

3 回答

  • 5

    由于您对系数施加的约束，您指定的模型不是线性模型，因此不能使用 lm 来拟合它 . 您需要使用非线性回归，例如 nls .

    > summary(nls(y ~ b0 + b1*x1 + b2*x2 + b1*b2*x1*x2, start=list(b0=0, b1=1, b2=1)))

Formula: y ~ b0 + b1 * x1 + b2 * x2 + b1 * b2 * x1 * x2

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b0 0.987203   0.049713   19.86   <2e-16 ***
b1 0.494438   0.007803   63.37   <2e-16 ***
b2 0.202396   0.003359   60.25   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1121 on 47 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 2.545e-06

    当您重新编写模型时，您确实可以看到模型是非线性的

    > summary(nls(y ~ b0+(1+b1*x1)*(1+b2*x2)-1, start=list(b0=0, b1=1, b2=1)))

Formula: y ~ b0 + (1 + b1 * x1) * (1 + b2 * x2) - 1

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b0 0.987203   0.049713   19.86   <2e-16 ***
b1 0.494438   0.007803   63.37   <2e-16 ***
b2 0.202396   0.003359   60.25   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1121 on 47 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 2.25e-06
    回复于
  • 1

    Brian提供了一种方法来拟合您指定的约束模型，但如果您对无约束模型比拟约束模型更合适感兴趣，可以使用delta方法来检验该假设 .

    # Let's make some fake data where the constrained model is true
n <- 100
b0 <- 2
b1 <- .2
b2 <- -1.3
b3 <- b1 * b2
sigma <- 1

x1 <- rnorm(n)
# make x1 and x2 correlated for giggles
x2 <- x1 + rnorm(n) 
# Generate data according to the model
y <- b0 + b1*x1 + b2*x2 + b3*x1*x2 + rnorm(n, 0, sigma)

# Fit full model y = b0 + b1*x1 + b2*x3 + b3*x1*x2 + error
o <- lm(y ~ x1 + x2 + x1:x2)

# If we want to do a hypothesis test of Ho: b3 = b1*b2
# this is the same as Ho: b3 - b1*b2 = 0
library(msm)
# Get estimate of the difference specified in the null
est <- unname(coef(o)["x1:x2"] - coef(o)["x1"] * coef(o)["x2"])
# Use the delta method to get a standard error for
# this difference
standerr <- deltamethod(~ x4 - x3*x2, coef(o), vcov(o))

# Calculate a test statistic.  We're relying on asymptotic
# arguments here so hopefully we have a decent sample size
z <- est/standerr
# Calculate p-value
pval <- 2 * pnorm(-abs(z))
pval

    我解释了delta方法的用途以及更多关于如何在this blog post中使用它的方法 .

    扩展Brian的答案，您可以通过将完整模型与约束模型进行比较来实现此目的 - 但是您必须使用nls来拟合整个模型才能轻松比较模型 .

    o2 <- nls(y ~ b0 + b1*x1 + b2*x2 + b1*b2*x1*x2, start=list(b0=0, b1=1, b2=1))
o3 <- nls(y ~ b0 + b1*x1 + b2*x2 + b3*x1*x2, start = list(b0 = 0, b1 = 1, b2 = 1, b3 = 1))
anova(o2, o3)
    回复于
  • 5

    有's no way to do what you'在 lm 要求，并且没有理由它能够做到这一点 . 您运行 lm 以获得系数的估计值 . 如果您没有在模型中包含预测变量 . 您可以使用 coef 提取所需的系数，然后将它们相乘 .

    请注意，离开交互是一个不同的模型，将产生不同的b1和b2 . 你也可以留下 I(x1 * x2) 而不使用系数 .

    至于你为什么要这样做，你的约束模型实际上比简单的加法模型更合适并不是一个好的先验理由 . 拥有更多的免费参数必然意味着模型更适合，但是你已经添加了一个约束，在现实世界中，它可以使它更适合 . 在这种情况下，你会认为它更好"baseline"用于与包括交互的模型进行比较吗？

    回复于

相关问题