Dopasowanie heteroscedastycznego uogólnionego modelu liniowego do odpowiedzi dwumianowych

9

Mam dane z następującego schematu eksperymentalnego: moje obserwacje są liczbami sukcesów ( K) z odpowiedniej liczby prób ( N), zmierzonymi dla dwóch grup, z których każda składa się z Iosobników, z Tzabiegów, w których w każdej takiej kombinacji czynników występują Rpowtórzenia . Stąd też zupełnie mam 2 * I * T * R K 'S i odpowiadającą ń ' s.

Dane pochodzą z biologii. Każdy osobnik jest genem, dla którego mierzę poziom ekspresji dwóch alternatywnych form (z powodu zjawiska zwanego alternatywnym składaniem). Zatem K jest poziomem ekspresji jednej z form, a N jest sumą poziomów ekspresji dwóch form. Zakłada się, że wybór między dwiema postaciami w jednym wyrażonym egzemplarzu jest eksperymentem Bernoulliego, stąd K z N.kopie następują po dwumianowym. Każda grupa składa się z ~ 20 różnych genów, a geny w każdej grupie mają pewną wspólną funkcję, która jest różna dla obu grup. Dla każdego genu w każdej grupie mam ~ 30 takich pomiarów z każdej z trzech różnych tkanek (zabiegów). Chcę oszacować wpływ, jaki grupa i leczenie mają na wariancję K / N.

Wiadomo, że ekspresja genów jest rozproszona, dlatego użycie ujemnego dwumianu w poniższym kodzie.

Np. RKod danych symulowanych:

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
}

Interesuje mnie oszacowanie wpływu, jaki grupa i leczenie mają na rozproszenie (lub wariancję) prawdopodobieństwa sukcesu (tj K/N.). Dlatego szukam odpowiedniego glm, w którym odpowiedź jest K / N, ale oprócz modelowania oczekiwanej wartości odpowiedzi modelowana jest również wariancja odpowiedzi.

Oczywiście na wariancję prawdopodobieństwa sukcesu dwumianowego ma wpływ liczba prób i leżące u ich podstaw prawdopodobieństwo sukcesu (im wyższa liczba prób i tym bardziej ekstremalne jest podstawowe prawdopodobieństwo sukcesu (tj. Bliskie 0 lub 1), tym niższy jest wariancja prawdopodobieństwa sukcesu), więc interesuje mnie głównie udział grupy i leczenia poza liczbą badań i leżącym u ich podstaw prawdopodobieństwem sukcesu. Myślę, że zastosowanie transformacji pierwiastka kwadratowego arcsin do odpowiedzi wyeliminuje tę drugą, ale nie liczbę prób.

Chociaż w symulowanym przykładzie dane powyżej projekt jest zrównoważony (równa liczba osobników w każdej z dwóch grup i identyczna liczba powtórzeń u każdej osoby z każdej grupy w każdym zabiegu), w moich rzeczywistych danych tak nie jest - obie grupy robią nie mają równej liczby osobników, a liczba powtórzeń jest różna. Wyobrażam sobie również, że jednostka powinna być ustawiona jako efekt losowy.

Wykreślenie wariancji próbki względem średniej próbki szacowanego prawdopodobieństwa sukcesu (oznaczonego jako p hat = K / N) każdego osobnika pokazuje, że skrajne prawdopodobieństwo sukcesu ma mniejszą wariancję: wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Eliminuje się to, gdy szacowane prawdopodobieństwo sukcesu transformowane jest za pomocą transformacji stabilizującej wariancję pierwiastka kwadratowego arcsin (oznaczonej jako arcsin (sqrt (p hat)): wprowadź opis zdjęcia tutaj

Wykreślenie wariancji próbki przekształconych oszacowanych prawdopodobieństw sukcesu w porównaniu ze średnią N pokazuje oczekiwaną ujemną zależność: wprowadź opis zdjęcia tutaj

Wykreślenie przykładowej wariancji przekształconych oszacowanych prawdopodobieństw sukcesu dla dwóch grup pokazuje, że grupa b ma nieco wyższe wariancje, i tak symulowałem dane: wprowadź opis zdjęcia tutaj

Na koniec, wykreślenie wariancji próbki przekształconych szacowanych prawdopodobieństw sukcesu dla trzech zabiegów nie pokazuje żadnej różnicy między zabiegami, tak symulowałem dane: wprowadź opis zdjęcia tutaj

Czy jest jakaś forma uogólnionego modelu liniowego, za pomocą którego mogę kwantyfikować wpływ grupy i leczenia na wariancję prawdopodobieństwa sukcesu?

Być może heteroscedastyczny uogólniony model liniowy lub jakaś forma logicznego modelu wariancji?

Coś w liniach modelu, który modeluje wariancję (y) = Zλ oprócz E (y) = Xβ, gdzie Z i X są odpowiednio regresorami odpowiednio średniej i wariancji, które w moim przypadku będą identyczne i obejmują leczenie (poziomy t.1, t.2 i t.3) i grupa (poziomy aib), i prawdopodobnie N i R, a zatem λ i β oszacują ich odpowiednie efekty.

Alternatywnie, mógłbym dopasować model do wariancji próbki między powtórzeniami każdego genu z każdej grupy w każdym traktowaniu, używając glm, który modeluje tylko oczekiwaną wartość odpowiedzi. Jedyne pytanie tutaj polega na tym, jak wyjaśnić fakt, że różne geny mają różną liczbę powtórzeń. Myślę, że wagi w glm mogłyby to wyjaśnić (wariancje próbek oparte na większej liczbie powtórzeń powinny mieć większą wagę), ale dokładnie jakie wagi należy ustawić?

Uwaga: próbowałem użyć dglmpakietu 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 

Efekt grupowy według dglm.fit jest dość słaby. Zastanawiam się, czy model jest dobrze ustawiony, czy też ma moc, jaką ma ten model.

użytkownik1701545
źródło
1
To jest trochę nieprzejrzyste. Czy możesz powiedzieć bardziej konkretnie, jaka jest Twoja sytuacja, być może konkretnie?
gung - Przywróć Monikę
1
Zredagowałem swój post. Ale ogólnie szukam glm, który oprócz E (y) = Xβ, modele Wariancja (y) = Zλ, gdzie kolumny X są regresorami dla średniej odpowiedzi, a kolumny Z są regresory wariancji odpowiedzi oraz β i λ są parametrami, które szacują ich wpływ odpowiednio na średnią i wariancję.
user1701545
2
Tyle zbieram, ale trudno jest to przemyśleć w tak mglistych kategoriach. Nie omawiasz żadnych zmiennych X ani Z, np. O czym jest twoje badanie? Jakie są twoje dane Itd.
gung - Przywróć Monikę
Nie ma problemu - zobacz zaktualizowany post.
user1701545

Odpowiedzi: