1

我知道这是一个众所周知的经典问题,但尽管我进行了研究,但我未能用我的数据解决这个问题。

我有这样的数据:

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 bloggerStackoverflow中的这些问题,我了解如何在添加-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 ?

任何有关以适当方式进行编程的建议都将不胜感激。

4

0 回答 0