如何比较R中的斜率答案

【问题标题】：How to compare slopes in R如何比较R中的斜率
【发布时间】：2017-12-01 10:57:25
【问题描述】：

我正在执行 ANCOVA，以测试每个性别（总变量，sexo）中体型（协变量，logLCC）和不同头部测量值（响应变量，logLP）之间的关系。我在lm 中得到了每个性别的斜率，我想将它们与 1 进行比较。更具体地说，我想知道斜率是显着高于还是小于 1，或者它们是否等于 1，如这将在它们的异速生长关系中具有不同的生物学意义。这是我的代码：

#Modelling my lm#
> lm.logLP.sexo.adu<-lm(logLP~logLCC*sexo, data=ADU)
> anova(lm.logLP.sexo.adu)
Analysis of Variance Table

Response: logLP
         Df Sum Sq Mean Sq  F value    Pr(>F)    
logLCC        1 3.8727  3.8727 3407.208 < 2.2e-16 ***
sexo          1 0.6926  0.6926  609.386 < 2.2e-16 ***
logLCC:sexo   1 0.0396  0.0396   34.829 7.563e-09 ***
Residuals   409 0.4649  0.0011                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Obtaining slopes#
> lm.logLP.sexo.adu$coefficients
 (Intercept)       logLCC        sexoM logLCC:sexoM 
  -0.1008891    0.6725818   -1.0058962    0.2633595 
> lm.logLP.sexo.adu1<-lstrends(lm.logLP.sexo.adu,"sexo",var="logLCC")
> lm.logLP.sexo.adu1
 sexo logLCC.trend         SE  df  lower.CL  upper.CL
 H       0.6725818 0.03020017 409 0.6132149 0.7319487
 M       0.9359413 0.03285353 409 0.8713585 1.0005241

Confidence level used: 0.95

#Comparing slopes#
> pairs(lm.logLP.sexo.adu1)
 contrast   estimate         SE  df t.ratio p.value
 H - M    -0.2633595 0.04462515 409  -5.902  <.0001

#Checking whether the slopes are different than 1#
#Computes Summary with statistics      
> s1<-summary(lm.logLP.sexo.adu)
> s1

Call:
lm(formula = logLP ~ logLCC * sexo, data = ADU)

Residuals:
 Min       1Q   Median       3Q      Max 
-0.13728 -0.02202 -0.00109  0.01880  0.12468 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.10089    0.12497  -0.807     0.42    
logLCC        0.67258    0.03020  22.271  < 2e-16 ***
sexoM        -1.00590    0.18700  -5.379 1.26e-07 ***
logLCC:sexoM  0.26336    0.04463   5.902 7.56e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03371 on 409 degrees of freedom
Multiple R-squared:  0.9083,    Adjusted R-squared:  0.9076 
F-statistic:  1350 on 3 and 409 DF,  p-value: < 2.2e-16

#Computes t-student H0: intercept=1. The estimation of coefficients and   their s.d. are in s1$coefficients
> t1<-(1-s1$coefficients[2,1])/s1$coefficients[2,2]
#Calculates two tailed probability
> pval<- 2 * pt(abs(t1), df = df.residual(lm.logLP.sexo.adu), lower.tail = FALSE)
> print(pval)
[1] 3.037231e-24

我在这里的几个线程中看到了整个过程。但我能理解的是，我的斜率与 1 只是不同。我如何检查它们是大于还是小于比1？

已编辑

解决了！

#performs one-side test H0=slope bigger than 1
pval<-pt(t1, df = df.residual(lm.logLP.sexo.adu), lower.tail = FALSE)
#performs one-side test H0=slope smaller than 1
pval<-pt(t1, df = df.residual(lm.logLP.sexo.adu), lower.tail = TRUE)

此外，应在单性别模型中进行测试。

【问题讨论】：

你在哪里看到你的斜率不同于 1？我可以看到比较两个斜率的测试的 p 值，这表明两个斜率不同。到目前为止没有与1的比较。我还可以看到 M 的斜率是 0.9359... 以及包含 1 [0.8713585, 1.0005241] 的置信区间。
在您提到的问题中，只有一个自变量，因此只有一个系数。你有 3 个系数，因为你有 logLCC、sex 和交互，所以你必须分别检查它们中的每一个。到目前为止，您所做的只是对您的第二个系数 (logLCC) 的测试。对于其余的系数，您必须这样做。(1-s1$coefficients[3,1])/s1$coefficients[3,2] 用于性系数等。
您正在进行的测试表明，您的所有系数在统计上都显着不同于 1。然后您可以查看系数的实际值并报告它们是否大于或小于 1。始终保持注意你在比较什么。例如，如果您的系数是 0.25，并且您发现它与 1 不同，则它小于 1。它不能更大:-)
另外，检查 2 tailed 和 1 tailed t 检验之间的差异。您正在执行的 2 尾检验是检验 coeff = 1。小的 p 值表明您有证据拒绝这个假设。因此，coeff 与 1 不同。然后检查实际的 coeff 值并报告它是更大还是更小。
您上面提到的系数来自函数lstrends，而不是来自模型。该模型不会为每种性别提供一个斜率，因为两种性别都属于一个变量（性别）。所以，模型总是会给你一个斜率。您可以为每种性别建立一个模型。您可以将性变量拆分为两个变量并运行另一个模型。你真的需要交互变量吗？我不熟悉该应用程序的研究。我建议您使用相同的应用程序遵循另一个类似的示例。

标签： r regression lm anova

【解决方案1】：

如何检查它们是大于还是小于 1？

如this post、this post 和您所讨论的那样，您可以进行 Wald 测试，您可以通过这些测试来计算

t1<-(1-s1$coefficients[2,1])/s1$coefficients[2,2]

或者，使用vcov 和coef 函数使代码更具可读性

fit <- lm.logLP.sexo.adu
t1<-(1-coef(fit)[1])/vcov(fit)[1, 1]

Wald 检验为您提供 t 统计量，可用于制作双面或单面 test。因此，您可以删除 abs 并根据要测试的尾部设置 lower.tail 参数。

【讨论】：