«Контрасты можно применять только к факторам, имеющим два и более уровней» — но все ли факторы имеют два и более уровней?
Раньше я успешно использовал пакет ClusterBootstrap в R для выполнения иерархической загрузки. Сейчас я пытаюсь применить его к новому набору данных и получаю ошибку при использовании функции clusbootglm для создания glm с использованием предикторов Pain и Drug на P7_Weight, которые используют Dam_Number (идентификатор помета) для контроля помета.
Вот dput для воссоздания рассматриваемого набора данных:
library(tidyverse)
library(ClusterBootstrap)
structure(list(Dam_Number = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L), levels =
c("6861",
"6867", "6868", "6870", "7273", "7274", "7275", "7276", "7300",
"7303", "7304"), class = "factor"), Pain_Drug = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels =
c("Sham_Saline",
"CCI_Oxy", "CCI_Saline", "Sham_Oxy"), class = "factor"), Pain =
structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), levels = c("Sham",
"CCI"), class = "factor"), Drug = structure(c(1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = c("Saline", "Oxy"),
class = "factor"),
Sex = c("M", "M", "M", "F", "F", "F", "F", "F", "F", "F",
"F", "F", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M",
"F", "F", "F", "F", "F", "F", "M", "M", "M", "M", "M", "M",
"F", "F", "F", "F", "M", "M", "M", "M", "M", "M", "F", "F",
"F", "F", "M", "M", "M", "M", "M", "F", "F", "F", "F", "F",
"F", "F", "M", "M", "M", "M", "M", "M", "M", "F", "F", "F",
"F", "F", "F", "F", "F", "M", "M", "M", "M", "M", "M", "M",
"M", "M", "F", "F", "F", "F", "F", "F", "F", "M", "M", "M",
"M", "M", "M", "M", "F", "F", "F", "F", "F", "M", "M", "M",
"M", "M", "M", "M", "F", "F", "F", "F", "F", "M", "F", "F",
"F", "F", "F", "F", "F", "M", "M", "M", "M", "M", "M", "M",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F"), Pup_Number =
c(1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17), P7_Weight
= c(14.2,
14.1, 14.3, 14, 13, 14.1, 14.1, 12.9, 13.6, 12.4, 14, 11.5,
12, 11, 11.6, 11.4, 8.7, 13.5, 10, 11.7, 10.9, 11.2, 11.3,
8.9, 11, 10.1, 10.8, 10, 16.5, 16.6, 14.5, 16.3, 14.9, 17.2,
15.8, 16.2, 15.7, 15.2, 15.4, 15.7, 16.3, 16.2, 16.3, 15.2,
16.4, 14.8, 15.2, 16.1, 14.8, 15.3, 14.6, 16, 15.9, 13.2,
15.2, 14.8, 14.8, 14.4, 13.2, 15.4, 12.3, 10.6, 11.5, 10,
12.8, 12.2, 12.7, 10.8, 10.8, 9.2, 10, 10, 10, 8.5, 10, 13.9,
16.5, 13, 16.3, 16.6, 13.8, 12.4, 14.3, 13.5, 15.2, 14.6,
13.1, 14, 15.8, 13.6, 11.7, 16, 14.1, 16.1, 14.2, 15.6, 16.3,
16, 14.5, 15.3, 15.5, 14.7, 13.7, 14.6, 14.1, 15.3, 14.6,
14.5, 13.1, 13.1, 14.2, 13.4, 12.1, 13.8, 13.6, 15.1, 13.2,
15, 15.7, 14.4, 14.7, 13.7, 13.5, 11.7, 12.3, 13.3, 10, 12.9,
12.9, 10, 10, 12.4, 11.3, 11.4, 13.2, 11.2, 12.7, 11.1, 9.5,
10), Ref = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), row.names = c(NA, -140L
), class = c("tbl_df", "tbl", "data.frame"))
Когда я запускаю функцию clusbootglm, как показано ниже:
model_weight <- clusbootglm(P7_Weight ~ Pain*Drug, data = weights, clusterid = Dam_Number, family = gaussian, B = 5000, confint.level = 0.95, n.cores = 1)
Он возвращает следующую ошибку (с включенной трассировкой):
'Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
11.stop("contrasts can be applied only to factors with 2 or more levels")
10.`contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]])
9.model.matrix.default(mt, mf, contrasts)
8.model.matrix(mt, mf, contrasts)
7.glm(model, family = family, data = data[obs, ])
6.coef(glm(model, family = family, data = data[obs, ]))
5.doTryCatch(return(expr), name, parentenv, handler)
4.tryCatchOne(expr, names, parentenv, handlers[[1L]])
3.tryCatchList(expr, classes, parentenv, handlers)
2.tryCatch(coef(glm(model, family = family, data = data[obs, ])),warning = function(x) rep(as.numeric(NA), p))
1.clusbootglm(P7_Weight ~ Pain * Drug, data = weights, clusterid = Dam_Number, family = gaussian, B = 5000, confint.level = 0.95, n.cores = 1)'
На сайте много дискуссий по поводу этой ошибки, но все вопросы сводятся к идее, что в наборе данных есть NA или постоянные переменные. В этом посте была отличная функция, позволяющая попытаться отладить эту ошибку. Мой набор данных с честью проходит эту отладку, заявляя, что все предикторы имеют более двух уровней:
'$nlevels
Dam_Number Pain_Drug Pain Drug Sex
11 4 2 2 2
$levels
$levels$Dam_Number
[1] "6861" "6867" "6868" "6870" "7273" "7274" "7275" "7276" "7300" "7303" "7304"
$levels$Pain_Drug
[1] "Sham_Saline" "CCI_Oxy" "CCI_Saline" "Sham_Oxy"
$levels$Pain
[1] "Sham" "CCI"
$levels$Drug
[1] "Saline" "Oxy"
$levels$Sex
[1] "F" "M"'
Интересно, что если я использую комбинированную переменную (Pain_Drug), имеющую 4 уровня, код работает без проблем.
model_weight <- clusbootglm(P7_Weight ~ Pain_Drug, data = weights, clusterid = Dam_Number, family = gaussian, B = 5000, confint.level = 0.95, n.cores = 1)
Кажется, есть проблема с аспектом взаимодействия.
Обычно это выглядит как другая ошибка, где указано количество образцов начальной загрузки NA. Это немного меньшая выборка, чем я использовал ранее, так что, возможно, это и так, но если я попытаюсь запустить модель с ЛИБО болью или лекарством вместо взаимодействия или комбинированной переменной, она все равно вернет ту же ошибку.
Итак, ваша модель основных эффектов работает Pain+Drug, но не модель взаимодействия Pain*Drug?
При начальной загрузке вы выбираете случайные переменные, и это может привести к пропуску некоторых уровней. Представьте, что у вас есть два уровня и все выборочные значения относятся к одной категории. Тогда вы получите эту ошибку (со мной это случалось не раз).
Если вы используете options(error = recover), то R остановится в той точке, где возникла ошибка. Когда я это сделал, это произошло в образце начальной загрузки 162. Предположительно, вы получите другую точку остановки с другим случайным начальным числом. В этот момент вы можете посмотреть переменные в вызове clusbootglm. (У меня есть n.cores == 1, что облегчает отладку.) Посмотрите data[obs, ].
Мне кажется, автор ClusterBootstrap упустил возможность возникновения этой ошибки. Это было бы легко исправить код: он уже ловит предупреждения, просто позвольте ему также ловить ошибки.
Ответы 1
Функция clusbootglm обрабатывает числовой ввод:
library(tidyverse)
library(ClusterBootstrap)
weights <- weights |>
mutate(Sex = as.factor(Sex)) |>
mutate(across(-Dam_Number, ~as.numeric(.))) |>
select(Dam_Number, Pain, Drug, P7_Weight)
model_weight <- clusbootglm(P7_Weight ~ Pain*Drug, data = weights,
clusterid = Dam_Number, family = gaussian, B = 5000, confint.level = 0.95, n.cores = 1)
Вероятно, это эквивалентная модель, но как сравниваются коэффициенты? Если я установлю B = 100, обе модели будут успешными, но параметризации совершенно разные.
Кажется, это решило проблему! Не знаю, почему сейчас требуются числовые значения, поскольку я вполне уверен, что использовал факторы в прошлом... Но это работает! Спасибо,
@HannahHarder: он работает, но дает ответы, отличные от тех, которые запрашивал исходный код. Знаете ли вы, как преобразовать одно в другое?
Я думаю, проблема в том, что в ваших кластерах нет никаких различий в факторах, поэтому в некоторых бутстрап-образцах один из факторов вообще не имеет никаких изменений. Я не знаю эту процедуру достаточно хорошо, чтобы предложить способ обойти эту проблему.