Ajustement d'un modèle linéaire généralisé hétéroscédastique pour les réponses binomiales

9

J'ai des données de la conception expérimentale suivante: mes observations sont des décomptes du nombre de succès ( K) sur le nombre correspondant d'essais ( N), mesurés pour deux groupes comprenant chacun des Iindividus, à partir de Ttraitements, où dans chaque combinaison de facteurs il y a des Rrépétitions . Par conséquent, au total, j'ai 2 * I * T * R K et les N correspondants .

Les données proviennent de la biologie. Chaque individu est un gène pour lequel je mesure le niveau d'expression de deux formes alternatives (en raison d'un phénomène appelé épissage alternatif). Par conséquent, K est le niveau d'expression d'une des formes et N est la somme des niveaux d'expression des deux formes. Le choix entre les deux formes dans une seule copie exprimée est supposé être une expérience de Bernoulli, donc K sur Ncopies suit un binôme. Chaque groupe est composé d'environ 20 gènes différents et les gènes de chaque groupe ont une fonction commune, qui est différente entre les deux groupes. Pour chaque gène dans chaque groupe, j'ai ~ 30 mesures de ce type dans chacun des trois tissus différents (traitements). Je veux estimer l'effet que le groupe et le traitement ont sur la variance de K / N.

L'expression des gènes est connue pour être sur-dispersée, d'où l'utilisation de binômes négatifs dans le code ci-dessous.

Par exemple, Rcode de données simulées:

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
}

Je souhaite estimer les effets du groupe et du traitement sur la dispersion (ou la variance) des probabilités de succès (c.-à K/N-d.). Par conséquent, je recherche un glm approprié dans lequel la réponse est K / N mais en plus de modéliser la valeur attendue de la réponse, la variance de la réponse est également modélisée.

De toute évidence, la variance d'une probabilité de succès binomiale est affectée par le nombre d'essais et la probabilité de succès sous-jacente (plus le nombre d'essais est élevé et plus la probabilité de succès sous-jacente est extrême (c.-à-d. Près de 0 ou 1), plus la variance de la probabilité de succès), je suis donc principalement intéressé par la contribution du groupe et du traitement au-delà de celle du nombre d'essais et de la probabilité de réussite sous-jacente. Je suppose que l'application de la transformation de la racine carrée de l'arcsin à la réponse éliminera cette dernière mais pas celle du nombre d'essais.

Bien que dans l'exemple de données simulées ci-dessus, la conception soit équilibrée (nombre égal d'individus dans chacun des deux groupes et nombre identique de répétitions dans chaque individu de chaque groupe dans chaque traitement), dans mes données réelles, ce n'est pas le cas - les deux groupes le font n'ont pas un nombre égal d'individus et le nombre de répétitions varie. En outre, j'imagine que l'individu devrait être défini comme un effet aléatoire.

Le tracé de la variance de l'échantillon par rapport à la moyenne de l'échantillon de la probabilité de réussite estimée (notée p hat = K / N) de chaque individu montre que les probabilités de succès extrêmes ont une variance plus faible: entrez la description de l'image ici entrez la description de l'image ici

Ceci est éliminé lorsque les probabilités de succès estimées sont transformées à l'aide de la transformation de stabilisation de la variance de la racine carrée d'arcsin (notée arcsin (sqrt (p hat)): entrez la description de l'image ici

Le tracé de la variance de l'échantillon des probabilités de succès estimées transformées par rapport à la moyenne N montre la relation négative attendue: entrez la description de l'image ici

Le tracé de la variance de l'échantillon des probabilités de succès estimées transformées pour les deux groupes montre que le groupe b a des variances légèrement plus élevées, c'est ainsi que j'ai simulé les données: entrez la description de l'image ici

Enfin, le traçage de la variance de l'échantillon des probabilités de succès estimées transformées pour les trois traitements ne montre aucune différence entre les traitements, c'est ainsi que j'ai simulé les données: entrez la description de l'image ici

Existe-t-il une forme de modèle linéaire généralisé permettant de quantifier les effets de groupe et de traitement sur la variance des probabilités de succès?

Peut-être un modèle linéaire généralisé hétéroscédastique ou une forme de modèle de variance log-linéaire?

Quelque chose dans les lignes d'un modèle qui modélise la variance (y) = Zλ en plus de E (y) = Xβ, où Z et X sont les régresseurs de la moyenne et de la variance, respectivement, qui dans mon cas seront identiques et incluront le traitement (niveaux t.1, t.2 et t.3) et le groupe (niveaux a et b), et probablement N et R, et donc λ et β estimeront leurs effets respectifs.

Alternativement, je pourrais adapter un modèle aux variances de l'échantillon entre les répliques de chaque gène de chaque groupe dans chaque traitement, en utilisant un glm qui ne modélise que la valeur attendue de la réponse. La seule question ici est de savoir comment tenir compte du fait que différents gènes ont différents nombres de répliques. Je pense que les poids dans un GLM pourraient expliquer cela (les variances d'échantillon qui sont basées sur plus de répétitions devraient avoir un poids plus élevé) mais exactement quels poids devraient être fixés?

Remarque: j'ai essayé d'utiliser le dglmpackage 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 

L'effet de groupe selon dglm.fit est assez faible. Je me demande si le modèle est bien réglé ou est la puissance de ce modèle.

user1701545
la source
1
C'est un peu opaque. Pouvez-vous nous dire plus précisément quelle est votre situation, peut-être concrètement?
gung - Rétablir Monica
1
J'ai édité mon message. Mais en général, je cherche un glm qui en plus de E (y) = Xβ, modélise la variance (y) = Zλ, où les colonnes de X sont les régresseurs de la moyenne de la réponse et les colonnes de Z sont les régresseurs de la variance de la réponse, et β et λ sont les paramètres qui estiment leurs effets sur la moyenne et la variance, respectivement.
user1701545
2
J'en comprends beaucoup, mais il est difficile de penser cela en des termes aussi nébuleux. Vous ne discutez d'aucune variable X ou Z, par exemple. Sur quoi porte votre étude? Quelles sont vos données? Etc.
gung - Réintégrer Monica
Pas de problème - voir l'article mis à jour.
user1701545

Réponses: