Anpassung eines heteroskedastischen verallgemeinerten linearen Modells für Binomialantworten

9

Ich habe Daten aus dem folgende experimentellen Design: meine Beobachtungen zählen die Anzahl der Erfolge sind ( K) aus Anzahl der Versuche von entsprechenden ( N), gemessen für zwei Gruppen von jeweils umfassten IIndividuen, von TBehandlungen, wobei in jeder solchen Faktorkombination gibt es RReplikate . Insgesamt habe ich also 2 * I * T * R K 's und entsprechende N ' s.

Die Daten stammen aus der Biologie. Jedes Individuum ist ein Gen, für das ich das Expressionsniveau von zwei alternativen Formen messe (aufgrund eines Phänomens, das als alternatives Spleißen bezeichnet wird). Daher ist K das Expressionsniveau einer der Formen und N ist die Summe der Expressionsniveaus der beiden Formen. Die Wahl zwischen den beiden Formen in einer einzelnen ausgedrückten Kopie wird als Bernoulli-Experiment angenommen, daher K aus N.Kopien folgen einem Binomial. Jede Gruppe besteht aus ~ 20 verschiedenen Genen und die Gene in jeder Gruppe haben eine gemeinsame Funktion, die sich zwischen den beiden Gruppen unterscheidet. Für jedes Gen in jeder Gruppe habe ich ~ 30 solcher Messungen von jedem von drei verschiedenen Geweben (Behandlungen). Ich möchte den Effekt abschätzen, den Gruppe und Behandlung auf die Varianz von K / N haben.

Es ist bekannt, dass die Genexpression überdispers ist, daher wird im folgenden Code ein negatives Binom verwendet.

ZB RCode simulierter Daten:

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
}

Ich bin daran interessiert, die Auswirkungen von Gruppe und Behandlung auf die Streuung (oder Varianz) der Erfolgswahrscheinlichkeiten (dh K/N) abzuschätzen . Daher suche ich nach einem geeigneten glm, in dem die Antwort K / N ist, aber zusätzlich zur Modellierung des erwarteten Werts der Antwort wird auch die Varianz der Antwort modelliert.

Die Varianz einer binomialen Erfolgswahrscheinlichkeit wird eindeutig von der Anzahl der Versuche und der zugrunde liegenden Erfolgswahrscheinlichkeit beeinflusst (je höher die Anzahl der Versuche ist und je extremer die zugrunde liegende Erfolgswahrscheinlichkeit ist (dh nahe 0 oder 1), desto geringer ist die Varianz der Erfolgswahrscheinlichkeit), daher interessiert mich hauptsächlich der Beitrag von Gruppe und Behandlung, der über die Anzahl der Studien und die zugrunde liegende Erfolgswahrscheinlichkeit hinausgeht. Ich denke, die Anwendung der Arcsin-Quadratwurzel-Transformation auf die Antwort wird letztere eliminieren, aber nicht die der Anzahl der Versuche.

Obwohl in den simulierten Beispieldaten oben das Design ausgewogen ist (gleiche Anzahl von Individuen in jeder der beiden Gruppen und identische Anzahl von Wiederholungen in jedem Individuum aus jeder Gruppe in jeder Behandlung), ist dies in meinen realen Daten nicht der Fall - die beiden Gruppen tun dies nicht die gleiche Anzahl von Personen haben und die Anzahl der Wiederholungen variiert. Ich würde mir auch vorstellen, dass die Person als zufälliger Effekt eingestellt werden sollte.

Das Auftragen der Stichprobenvarianz gegen den Stichprobenmittelwert der geschätzten Erfolgswahrscheinlichkeit (bezeichnet als p hat = K / N) jedes Individuums zeigt, dass extreme Erfolgswahrscheinlichkeiten eine geringere Varianz aufweisen: Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein

Dies wird beseitigt, wenn die geschätzten Erfolgswahrscheinlichkeiten unter Verwendung der Arcsin-Quadratwurzel-Varianz-Stabilisierungstransformation (bezeichnet als Arcsin (sqrt (p hat)) transformiert werden: Geben Sie hier die Bildbeschreibung ein

Das Auftragen der Stichprobenvarianz der transformierten geschätzten Erfolgswahrscheinlichkeiten gegen den Mittelwert N zeigt die erwartete negative Beziehung: Geben Sie hier die Bildbeschreibung ein

Das Auftragen der Stichprobenvarianz der transformierten geschätzten Erfolgswahrscheinlichkeiten für die beiden Gruppen zeigt, dass Gruppe b geringfügig höhere Varianzen aufweist. So habe ich die Daten simuliert: Geben Sie hier die Bildbeschreibung ein

Schließlich zeigt das Auftragen der Stichprobenvarianz der transformierten geschätzten Erfolgswahrscheinlichkeiten für die drei Behandlungen keinen Unterschied zwischen den Behandlungen. So habe ich die Daten simuliert: Geben Sie hier die Bildbeschreibung ein

Gibt es irgendeine Form eines verallgemeinerten linearen Modells, mit dem ich die Gruppen- und Behandlungseffekte auf die Varianz der Erfolgswahrscheinlichkeiten quantifizieren kann?

Vielleicht ein heteroskedastisches verallgemeinertes lineares Modell oder irgendeine Form eines loglinearen Varianzmodells?

Etwas in den Linien eines Modells, das zusätzlich zu E (y) = Xβ die Varianz (y) = Zλ modelliert, wobei Z und X die Regressoren des Mittelwerts bzw. der Varianz sind, die in meinem Fall identisch sind und umfassen Behandlung (Niveaus t.1, t.2 und t.3) und Gruppe (Niveaus a und b) und wahrscheinlich N und R, und daher werden λ und β ihre jeweiligen Wirkungen abschätzen.

Alternativ könnte ich ein Modell an die Probenvarianzen über Replikate jedes Gens aus jeder Gruppe in jeder Behandlung anpassen, wobei ein glm verwendet wird, das nur den erwarteten Wert der Antwort modelliert. Die einzige Frage hier ist, wie die Tatsache zu erklären ist, dass verschiedene Gene unterschiedliche Anzahlen von Replikaten haben. Ich denke, die Gewichte in einem glm könnten das erklären (Stichprobenvarianzen, die auf mehr Replikaten basieren, sollten ein höheres Gewicht haben), aber welche Gewichte sollten genau eingestellt werden?

Hinweis: Ich habe versucht, das dglmR-Paket zu verwenden:

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 

Der Gruppeneffekt ist laut dglm.fit ziemlich schwach. Ich frage mich, ob das Modell richtig eingestellt ist oder welche Leistung dieses Modell hat.

user1701545
quelle
1
Das ist etwas undurchsichtig. Können Sie genauer sagen, wie Ihre Situation ist, vielleicht konkret?
Gung - Reinstate Monica
1
Ich habe meinen Beitrag bearbeitet. Aber im Allgemeinen suche ich nach einem glm, das zusätzlich zu E (y) = Xβ die Varianz (y) = Zλ modelliert, wobei die Spalten von X die Regressoren für den Mittelwert der Antwort und die Spalten von Z sind Die Regressoren für die Varianz der Antwort und β und λ sind die Parameter, die ihre Auswirkungen auf den Mittelwert bzw. die Varianz abschätzen.
user1701545
2
Ich sammle so viel, aber es ist schwierig, dies so nebulös zu überdenken. Sie diskutieren keine X- oder Z-Variablen, z. Worum geht es in Ihrer Studie? Was sind Ihre Daten? Usw.
gung - Reinstate Monica
Kein Problem - siehe aktualisierten Beitrag.
user1701545

Antworten: