1

我正在运行线性模型来查看所涉及的独立因素的重要性。示例模型是：`

mymod1 <- lm(temp ~ bgrp+psex+tb,data=mydat)
summary(mymod1)`

我查看摘要以检查每个因素的重要性：

lm(formula = temp ~ bgrp + psex + tb, data = mydat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6877 -0.2454  0.0768  0.3916  1.6561 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.324459   0.186081 200.581  < 2e-16 ***
bgrp         0.256794   0.066167   3.881 0.000115 ***
psex         0.144669   0.055140   2.624 0.008913 ** 
tb           0.019818   0.009342   2.121 0.034287 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6888 on 621 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.03675,    Adjusted R-squared: 0.03209 
F-statistic: 7.897 on 3 and 621 DF,  p-value: 3.551e-05

现在，我想看看bgrp（1和2）和psex（1和2）这两个级别的解决方案。

如果您能帮我解决这个问题，我将不胜感激。

提前谢谢你，

巴兹

编辑：

我运行了您建议的第一个模型并得到以下结果：

mydat$bgrp <- as.factor(mydat$bgrp)

> summary(lm(temp ~ bgrp+psex+tb-1,data=mydat))

Call:
lm(formula = temp ~ bgrp + psex + tb - 1, data = apirt)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6877 -0.2454  0.0768  0.3916  1.6561 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
bgrp1 37.725922   0.135486 278.449  < 2e-16 ***
bgrp2 37.982716   0.129558 293.171  < 2e-16 ***
psex2  0.144669   0.055140   2.624  0.00891 ** 
tb     0.019818   0.009342   2.121  0.03429 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6888 on 621 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.9997,     Adjusted R-squared: 0.9997 
F-statistic: 4.788e+05 on 4 and 621 DF,  p-value: < 2.2e-16

从上面的系数表来看，bgrp1 和 bgrp2 似乎是有道理的：bgrp1 代表产仔数较大的母系，后代较轻，这导致后代的直肠温度较低（37.70 摄氏度）。另一方面，bgrp2 代表产仔数较小、后代较重的终端线，这导致直肠温度较高（37.98 摄氏度）。我只是想知道，是否可以对 psex1 和 psex2 执行相同的操作，但是系数表中显示的内容可能是由于您之前所说的。

编辑：嗨，马克，

我尝试了您建议的两个选项，我可以看到 bgrp1 和 psex1 具有相同的值：

> mybgrp <- lm(formula = temp ~ bgrp+psex+tb-1, data = mydat)
> mybgrp

Call:
lm(formula = temp ~ bgrp + psex + tb - 1, data = mydat)

Coefficients:
   bgrp1     bgrp2     psex2        tb  
37.72592  37.98272   0.14467   0.01982  

> summary(mybgrp)

Call:
lm(formula = temp ~ bgrp + psex + tb - 1, data = mydat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6877 -0.2454  0.0768  0.3916  1.6561 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
bgrp1 37.725922   0.135486 278.449  < 2e-16 ***
bgrp2 37.982716   0.129558 293.171  < 2e-16 ***
psex2  0.144669   0.055140   2.624  0.00891 ** 
tb     0.019818   0.009342   2.121  0.03429 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6888 on 621 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.9997,     Adjusted R-squared: 0.9997 
F-statistic: 4.788e+05 on 4 and 621 DF,  p-value: < 2.2e-16 

> mypsex <- lm(formula = temp ~ psex+bgrp+tb-1, data = mydat)
> mypsex

Call:
lm(formula = temp ~ psex + bgrp + tb - 1, data = mydat)

Coefficients:
   psex1     psex2     bgrp2        tb  
37.72592  37.87059   0.25679   0.01982  

> summary(mypsex)

Call:
lm(formula = temp ~ psex + bgrp + tb - 1, data = mydat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6877 -0.2454  0.0768  0.3916  1.6561 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
psex1 37.725922   0.135486 278.449  < 2e-16 ***
psex2 37.870591   0.135908 278.649  < 2e-16 ***
bgrp2  0.256794   0.066167   3.881 0.000115 ***
tb     0.019818   0.009342   2.121 0.034287 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6888 on 621 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.9997,     Adjusted R-squared: 0.9997 
F-statistic: 4.788e+05 on 4 and 621 DF,  p-value: < 2.2e-16

谢谢！

4

1 回答 1

1

如果变量只有两个水平（1 vs 2），则与（0 vs 1）相同，斜率是这两个水平之一。变量的另一个级别包含在截距项中。

也许试试

lm(formula = temp ~ bgrp + psex + tb - 1 , data = mydat)

删除拦截，看看是否给你你想要的。

再说一次，也许我没有正确理解你的问题。

编辑：

当我使用假数据并设置

bgrp <- as.factor(bgrp)
psex <- as.factor(psex)

如果没有截距，我会得到 2 个因素之一的 2 个斜率。我相信 R 保持第二个因素 = 0 的第二个斜率。

编辑2：

该模型将为 bgrp 和 psex 的每种组合提供单独的斜率。该模型包括 bgrp 和 psex 之间的交互，然后删除截距以及 bgrp 和 psex 主效应：

mymod2 <- lm(temp ~ bgrp + psex + bgrp * psex + tb - 1 - bgrp - psex)

编辑3：

如果您习惯使用 SAS 并尝试在 SAS 和 R 中运行相同的分析，您可能会发现这两个程序最初不会返回相同的估计值。这可能是因为 SAS 和 R 默认为截距选择不同的因子水平。您可以更改 R 中截距的默认因子级别以与 SAS 使用的相匹配，然后您可能会发现两个程序给出相同的答案。

将以下 R 代码与此处的 SAS 输出进行比较：

http://support.sas.com/kb/38/384.html

SAS代码使用选项“解决方案”的地方：

my.data <- matrix(c(
'A', 'F',   9, 25,  
'A', 'F',   3, 19,  
'A', 'F',   4, 18,  
'A', 'F',  11, 28,  
'A', 'F',   7, 23,
'A', 'M',  11, 27,  
'A', 'M',   9, 24,  
'A', 'M',   9, 25,  
'A', 'M',  10, 28,  
'A', 'M',  10, 26,
'D', 'F',   4, 37,  
'D', 'F',  12, 54,  
'D', 'F',   3, 33,  
'D', 'F',   6, 41,  
'D', 'F',   9, 47,
'D', 'M',   5, 36,  
'D', 'M',   4, 36,  
'D', 'M',   7, 40,  
'D', 'M',  10, 46,  
'D', 'M',   8, 42,
'G', 'F',  10, 70,  
'G', 'F',  11, 75,  
'G', 'F',   7, 60,  
'G', 'F',   9, 69,  
'G', 'F',  10, 71,
'G', 'M',   3, 47,  
'G', 'M',   8, 60,  
'G', 'M',  11, 70,  
'G', 'M',   4, 49,  
'G', 'M',   4, 50
), nrow = 30, byrow=T, 
dimnames = list(NULL, c("drug","gender","x","y")));


my.data <- as.data.frame(my.data, stringsAsFactors=F)
my.data

my.data$y      <- as.numeric(my.data$y)
my.data$x      <- as.numeric(my.data$x)
my.data$drug   <- as.factor(my.data$drug)
my.data$gender <- as.factor(my.data$gender)

str(my.data)

my.data$drug   <- relevel(my.data$drug, ref="G")
my.data$gender <- relevel(my.data$gender, ref="M")



my.mod1 <- lm(my.data$y ~ my.data$drug)
my.mod1
summary(my.mod1)

my.mod2 <- lm(my.data$y ~ my.data$drug-1)
my.mod2
summary(my.mod2)

my.mod3 <- lm(my.data$y ~ my.data$drug + my.data$gender + 
                          my.data$drug * my.data$gender )
my.mod3
summary(my.mod3)

my.mod4 <- lm(my.data$y ~ my.data$drug + my.data$gender + 
                          my.data$drug * my.data$gender - 1 )
my.mod4
summary(my.mod4)

my.mod5 <- lm(my.data$y ~ my.data$drug + my.data$x + 
                          my.data$drug * my.data$x )
my.mod5
summary(my.mod5)

my.mod6 <- lm(my.data$y ~ my.data$drug + my.data$x + 
                          my.data$drug * my.data$x - 1 )
my.mod6
summary(my.mod6)
于 2012-03-11T19:46:11.953 回答