Подгонка гетероскедастической обобщенной линейной модели для биномиальных ответов

9

У меня есть данные из следующего экспериментального плана: мои наблюдения - это подсчет числа успехов ( K) из соответствующего числа испытаний ( N), измеренных для двух групп, каждая из которых состоит из Iиндивидуумов, из Tобработок, где в каждой такой комбинации факторов есть Rповторения , Таким образом, в целом у меня есть 2 * I * Т * Р К «ы и соответствующих N » с.

Данные из биологии. Каждый индивид - это ген, для которого я измеряю уровень экспрессии двух альтернативных форм (из-за явления, называемого альтернативным сплайсингом). Следовательно, K - уровень экспрессии одной из форм, а N - сумма уровней экспрессии двух форм. Предполагается, что выбор между двумя формами в одной выраженной копии является экспериментом Бернулли, следовательно, K из NКопии следует за бином. Каждая группа состоит из ~ 20 различных генов, и у генов в каждой группе есть некоторая общая функция, которая отличается между двумя группами. Для каждого гена в каждой группе у меня есть ~ 30 таких измерений для каждой из трех разных тканей (обработок). Я хочу оценить влияние группы и лечения на дисперсию K / N.

Известно, что экспрессия генов избыточна, поэтому в приведенном ниже коде используется отрицательный бином.

Например, Rкод моделируемых данных:

library(MASS)
set.seed(1)
I = 20 # individuals in each group
G = 2  # groups
T = 3  # treatments
R = 30 # replicates of each individual, in each group, in each treatment

groups     = letters[1:G]
ids        = c(sapply(groups, function(g){ paste(rep(g, I), 1:I, sep=".") }))
treatments = paste(rep("t", T), 1:T, sep=".")
 # create random mean number of trials for each individual and 
 #  dispersion values to simulate trials from a negative binomial:
mean.trials = rlnorm(length(ids), meanlog=10, sdlog=1)
thetas      = 10^6/mean.trials
 # create the underlying success probability for each individual:
p.vec = runif(length(ids), min=0, max=1)
 # create a dispersion factor for each success probability, where the 
 #  individuals of group 2 have higher dispersion thus creating a group effect:
dispersion.vec = c(runif(length(ids)/2, min=0, max=0.1),
                   runif(length(ids)/2, min=0, max=0.2))
 # create empty an data.frame:
data.df = data.frame(id=rep(sapply(ids, function(i){ rep(i, R) }), T),
                     group=rep(sapply(groups, function(g){ rep(g, I*R) }), T),
                     treatment=c(sapply(treatments, 
                                        function(t){ rep(t, length(ids)*R) })),
                     N=rep(NA, length(ids)*T*R), 
                     K=rep(NA, length(ids)*T*R) )
 # fill N's and K's - trials and successes
for(i in 1:length(ids)){
  N     = rnegbin(T*R, mu=mean.trials[i], theta=thetas[i])
  probs = runif(T*R, min=max((1-dispersion.vec[i])*p.vec[i],0),
                max=min((1+dispersion.vec)*p.vec[i],1))
  K     = rbinom(T*R, N, probs)
  data.df$N[which(as.character(data.df$id) == ids[i])] = N
  data.df$K[which(as.character(data.df$id) == ids[i])] = K
}

Мне интересно оценить влияние, которое группа и лечение оказывают на дисперсию (или дисперсию) вероятностей успеха (т. Е. K/N). Поэтому я ищу подходящий glm, в котором ответом является K / N, но в дополнение к моделированию ожидаемого значения ответа также моделируется дисперсия ответа.

Ясно, что на дисперсию вероятности биномиального успеха влияет количество испытаний и основная вероятность успеха (чем больше число испытаний и чем более экстремальной является основная вероятность успеха (т. Е. Около 0 или 1), тем ниже дисперсия вероятности успеха), поэтому я в основном интересуюсь вкладом группы и лечения помимо количества испытаний и основной вероятности успеха. Я предполагаю, что применение преобразования квадратного корня arcsin к ответу устранит последнее, но не количество испытаний.

Хотя в приведенном выше примере смоделированных данных дизайн сбалансирован (равное количество индивидуумов в каждой из двух групп и одинаковое количество повторов у каждого индивидуума из каждой группы в каждом лечении), в моих реальных данных это не так - две группы не имеют равное количество особей и количество повторов варьируется. Кроме того, я предполагаю, что индивидуум должен быть установлен как случайный эффект.

Построение выборочной дисперсии по отношению к выборочному среднему значению предполагаемой вероятности успеха (обозначается как p hat = K / N) каждого человека показывает, что экстремальные вероятности успеха имеют меньшую дисперсию: введите описание изображения здесь введите описание изображения здесь

Это устраняется, когда предполагаемые вероятности успеха преобразуются с использованием преобразования стабилизации дисперсии квадратного корня arcsin (обозначается как arcsin (sqrt (p hat)): введите описание изображения здесь

Построение выборочной дисперсии преобразованных оценочных вероятностей успеха в зависимости от среднего N показывает ожидаемые отрицательные отношения: введите описание изображения здесь

Построение выборочной дисперсии преобразованных оценочных вероятностей успеха для двух групп показывает, что группа b имеет несколько более высокие дисперсии, как я смоделировал данные: введите описание изображения здесь

Наконец, построение выборочной дисперсии преобразованных оценочных вероятностей успеха для трех обработок не показывает разницы между обработками, как я смоделировал данные: введите описание изображения здесь

Существует ли какая-либо форма обобщенной линейной модели, с помощью которой я могу количественно оценить влияние группы и лечения на дисперсию вероятностей успеха?

Возможно, гетероскедастическая обобщенная линейная модель или некоторая форма логлинейной модели дисперсии?

Что-то в строках модели, которая моделирует дисперсию (y) = Zλ в дополнение к E (y) = Xβ, где Z и X - регрессоры среднего значения и дисперсии соответственно, которые в моем случае будут идентичны и включают обработка (уровни t.1, t.2 и t.3) и группа (уровни a и b), и, вероятно, N и R, и, следовательно, λ и β оценят их соответствующие эффекты.

В качестве альтернативы, я мог бы подобрать модель для выборочных дисперсий между повторностями каждого гена из каждой группы в каждой обработке, используя glm, который только моделирует ожидаемое значение ответа. Единственный вопрос здесь заключается в том, как объяснить тот факт, что разные гены имеют разное количество повторностей. Я думаю, что веса в glm могли бы объяснить это (выборочные отклонения, которые основаны на большем количестве повторов, должны иметь больший вес), но какие именно веса должны быть установлены?

Примечание: я пробовал использовать dglmпакет R:

library(dglm)
dglm.fit = dglm(formula = K/N ~ 1, dformula = ~ group + treatment, family = quasibinomial, weights = N, data = data.df)
summary(dglm.fit)
Call: dglm(formula = K/N ~ 1, dformula = ~group + treatment, family = quasibinomial, 
    data = data.df, weights = N)

Mean Coefficients:
               Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -0.09735366 0.01648905 -5.904138 3.873478e-09
(Dispersion Parameters for quasibinomial family estimated as below )

    Scaled Null Deviance: 3600 on 3599 degrees of freedom
Scaled Residual Deviance: 3600 on 3599 degrees of freedom

Dispersion Coefficients:
                Estimate Std. Error      z value  Pr(>|z|)
(Intercept)  9.140517930 0.04409586 207.28746254 0.0000000
group       -0.071009599 0.04714045  -1.50634107 0.1319796
treatment   -0.001469108 0.02886751  -0.05089138 0.9594121
(Dispersion parameter for Gamma family taken to be 2 )

    Scaled Null Deviance: 3561.3 on 3599 degrees of freedom
Scaled Residual Deviance: 3559.028 on 3597 degrees of freedom

Minus Twice the Log-Likelihood: 29.44568 
Number of Alternating Iterations: 5 

Групповой эффект согласно dglm.fit довольно слабый. Интересно, установлена ​​ли модель правильно или это сила, которую имеет эта модель.

user1701545
источник
1
Это немного непрозрачно. Можете ли вы сказать более конкретно, какова ваша ситуация, возможно, в конкретных терминах?
gung - Восстановить Монику
1
Я отредактировал свой пост. Но в целом я ищу glm, который в дополнение к E (y) = Xβ, модели Variance (y) = Zλ, где столбцы X являются регрессорами для среднего ответа, а столбцы Z регрессоры для дисперсии ответа, а β и λ - параметры, которые оценивают их влияние на среднее значение и дисперсию соответственно.
user1701545
2
Я так много собираю, но об этом трудно думать в таких туманных выражениях. Вы не обсуждаете переменные X или Z, например. О чем твоя учеба? Каковы ваши данные? И т.д.
Ган - Восстановить Монику
Нет проблем - смотрите обновленный пост.
user1701545

Ответы: