我知道这是一个众所周知的经典问题,但尽管我进行了研究,但我未能用我的数据解决这个问题。
我有这样的数据:
df
SNP Site Color Frequence
1 scaffold10000|size69197_10061 K Green 0.4404348
2 scaffold10000|size69197_10061 G Green 0.6700000
3 scaffold10000|size69197_10061 G Red 0.7171429
4 scaffold10000|size69197_10061 K Yellow 0.7937500
5 scaffold10000|size69197_10061 T Yellow 0.7202174
6 scaffold10000|size69197_10061 E Red 0.7373469
7 scaffold10000|size69197_10061 G Yellow 0.6150000
8 scaffold10000|size69197_10061 T Red 0.5668750
9 scaffold10000|size69197_10061 K Red 0.6190385
10 scaffold10000|size69197_10061 T Green 0.5629412
11 scaffold10000|size69197_10061 E Yellow 0.8312500
12 scaffold10000|size69197_10061 E Green 0.5474286
我想知道这个 SNP(名为“scaffold10000|size69197_10061”)的三种颜色和四个位点之间是否存在统计差异。我想考虑这些变量(3 种颜色和 4 个站点),这就是我选择glm()
函数的原因。
model <- glm(formula = Frequence ~ Color + Site, family=quasibinomial(), data=df)
这给了我这些系数:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2905 0.3105 0.936 0.3856
ColorRed 0.4450 0.3084 1.443 0.1991
ColorYellow 0.8298 0.3215 2.581 0.0417 *
SiteT -0.2268 0.3644 -0.622 0.5566
SiteK -0.2221 0.3645 -0.609 0.5646
SiteE 0.1809 0.3760 0.481 0.6475
---
所以没有出现绿色和 G 站点(因为如果我理解正确,两者都是分类的)。
根据R blogger和Stackoverflow中的这些问题,我了解如何在添加-1
或+ 0
公式中删除截距(以使模型更易于理解)。
model <- glm(formula = Frequence ~ Color + Site - 1, family=quasibinomial(), data=df)
所以我至少出现了一个分类变量:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
ColorGreen 0.2905 0.3105 0.936 0.3856
ColorRed 0.7355 0.3185 2.309 0.0603 .
ColorYellow 1.1202 0.3319 3.376 0.0149 *
SiteT -0.2268 0.3644 -0.622 0.5566
SiteK -0.2221 0.3645 -0.609 0.5646
SiteE 0.1809 0.3760 0.481 0.6475
为了出现第 4 个站点,我未能对某些内容进行编码
首先,我尝试合并 2 个不同的模型:
model1 <- glm(formula = Frequence ~ Site - 1, family=quasibinomial(), data=df)
model2 <- glm(formula = Frequence ~ Color - 1, family=quasibinomial(), data=df)
以不同的方式但没有工作(并且可能没有意义..)
把其他人-1
或+ 0
都没有工作:
model <- glm(formula = Frequence ~ 0 + Color + Site - 1, family=quasibinomial(), data=df)
根据这个类似问题的答案(以及这个关于lm()
,只需在参数上添加两个总和为零的约束:
contrasts(ok$Site) <- contr.sum(4, contrasts=F)
contrasts(ok$Color) <- contr.sum(3, contrasts=F)
或使用这个(我不记得在 Stackoverflow 上的每一步)
relevel(ok$Site, "E")
relevel(ok$Site, "T")
relevel(ok$Site, "K")
relevel(ok$Site, "G")
并重新运行模型。但这两种可能性也都失败了。
所以我尝试拆分data.frame
以便手动将变量放入模型中:
df2
SNP Site Color Frequence Green Yellow Red K G T E
1 scaffold10000|size69197_10061 K Green 0.4404348 TRUE FALSE FALSE TRUE FALSE FALSE FALSE
2 scaffold10000|size69197_10061 G Green 0.6700000 TRUE FALSE FALSE FALSE TRUE FALSE FALSE
3 scaffold10000|size69197_10061 G Red 0.7171429 FALSE FALSE TRUE FALSE TRUE FALSE FALSE
4 scaffold10000|size69197_10061 K Yellow 0.7937500 FALSE TRUE FALSE TRUE FALSE FALSE FALSE
5 scaffold10000|size69197_10061 T Yellow 0.7202174 FALSE TRUE FALSE FALSE FALSE TRUE FALSE
6 scaffold10000|size69197_10061 E Red 0.7373469 FALSE FALSE TRUE FALSE FALSE FALSE TRUE
7 scaffold10000|size69197_10061 G Yellow 0.6150000 FALSE TRUE FALSE FALSE TRUE FALSE FALSE
8 scaffold10000|size69197_10061 T Red 0.5668750 FALSE FALSE TRUE FALSE FALSE TRUE FALSE
9 scaffold10000|size69197_10061 K Red 0.6190385 FALSE FALSE TRUE TRUE FALSE FALSE FALSE
10 scaffold10000|size69197_10061 T Green 0.5629412 TRUE FALSE FALSE FALSE FALSE TRUE FALSE
11 scaffold10000|size69197_10061 E Yellow 0.8312500 FALSE TRUE FALSE FALSE FALSE FALSE TRUE
12 scaffold10000|size69197_10061 E Green 0.5474286 TRUE FALSE FALSE FALSE FALSE FALSE TRUE
(TRUE 和 FALSE 可以用 . 变为 0 和 1 df2[df2=="FALSE"]<-0
。
model <- glm(formula=Frequence ~ Red + Green + Yellow + K + T + E + G -1, family=quasibinomial(), data=df2)
现在,所有变量都在系数中:
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
RedFALSE 1.1202 0.3319 3.376 0.0149 *
RedTRUE 0.7355 0.3185 2.309 0.0603 .
GreenTRUE -0.8298 0.3215 -2.581 0.0417 *
YellowTRUE NA NA NA NA
KTRUE -0.2221 0.3645 -0.609 0.5646
TTRUE -0.2268 0.3644 -0.622 0.5566
ETRUE 0.1809 0.3760 0.481 0.6475
GTRUE NA NA NA NA
但是NA
现在出现了。
根据Stackexchange中的这个问题,我检查了模型矩阵是否有满秩,答案是否定的。
# Get model matrix ...
X <- model.matrix(~ Red + Green + Yellow + K + T + E + G - 1, family=quasibinomial(), data=as.data.frame(ok))
> X
RedFALSE RedTRUE GreenTRUE YellowTRUE KTRUE TTRUE ETRUE GTRUE
1 1 0 1 0 1 0 0 0
2 1 0 1 0 0 0 0 1
3 0 1 0 0 0 0 0 1
4 1 0 0 1 1 0 0 0
5 1 0 0 1 0 1 0 0
6 0 1 0 0 0 0 1 0
7 1 0 0 1 0 0 0 1
8 0 1 0 0 0 1 0 0
9 0 1 0 0 1 0 0 0
10 1 0 1 0 0 1 0 0
11 1 0 0 1 0 0 1 0
12 1 0 1 0 0 0 1 0
# Get rank of model matrix
qr(X)$rank
> 6
# Get number of parameters of the model = number of columns of model matrix
ncol(X)
> 8
因此,如果没有-1
,则第一列X
是截距,如果有,-1
则红色列是重复的(一个代表 TRUE,一个代表 FALSE)。
因此,有 8 列和 6 个等级。通常我应该有 14 列和 14 列不?(7 个变量(3 个颜色和 4 个站点)* 2(真或假))
那么,如何对模型进行编码以强制获取所有变量的 Pvalues ?
任何有关以适当方式进行编程的建议都将不胜感激。