Я понимаю, что это хорошо известная классическая проблема, но, несмотря на мои исследования, мне не удалось исправить эту проблему с помощью моих данных.
У меня есть такие данные:
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
---
Так что сайт Green и 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
(ИСТИНА и ЛОЖЬ можно изменить на 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, красный столбец дублируется (один для ИСТИНЫ и один для ЛОЖИ).
Итак, есть 8 столбцов и 6 рангов.
Обычно у меня должно было быть 14 столбцов, а 14 - нет? (7 переменных (3 цвета и 4 сайта) * 2 (ИСТИНА или ЛОЖЬ))
Итак, как мне кодировать мою модель, чтобы принудительно получать значения Pvalue для всех переменных?
Любые советы по правильному программированию, они будут очень благодарны.