Tengo los datos de la siguiente diseño experimental: mis observaciones se cuenta el número de éxitos (K) de la correspondiente número de ensayos (N), medido por dos grupos, cada uno compuesto de I a los individuos, a partir de T tratamientos, donde en cada combinación de factor hay R repeticiones. Por lo tanto, en total tengo 2 * I * T * R K's y su correspondiente N's.
Los datos son de la biología. Cada individuo es un gen que puedo medir el nivel de expresión de las dos formas alternativas (debido a un fenómeno llamado splicing alternativo). Por lo tanto, K es el nivel de expresión de uno de los formularios y N es la suma de los niveles de expresión de las dos formas. La elección entre las dos formas en una sola expresó copia se supone que es un experimento de Bernoulli, por lo tanto K de N copias sigue una binomial. Cada grupo está compuesto de ~20 diferentes genes y los genes de cada grupo tienen en común la función, que es diferente entre los dos grupos. Para cada gen en cada grupo de he ~30 de tales mediciones de cada uno de los tres tejidos diferentes (tratamientos). Yo quiero estimar el efecto de que el grupo de tratamiento en la variación de K/N.
La expresión génica es conocido por ser overdispersed por lo tanto el uso de la binomial negativa en el código de abajo.
E. g., R código de datos simulados:
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
}
Estoy interesado en la estimación de los efectos de ese grupo y de trato en la dispersión (o la varianza) de las probabilidades de éxito (es decir, K/N). Por lo tanto, estoy en busca de un adecuado glm en el que la respuesta es K/N pero además de modelar el valor esperado de la respuesta de la varianza de la respuesta también es modelo.
Claramente, la varianza de una binomial probabilidad de éxito es afectado por el número de ensayos subyacente y la probabilidad de éxito (el más alto es el número de ensayos y la más extrema de las subyacentes de la probabilidad de éxito es (es decir, cerca de 0 o 1), menor es la variación de la probabilidad de éxito), por lo que estoy principalmente interesado en la contribución del grupo de tratamiento y más allá de que el número de ensayos subyacente y la probabilidad de éxito. Supongo que la aplicación de la arcsen transformación de raíz cuadrada a la respuesta, se eliminará el último pero no el de el número de ensayos.
A pesar de que en la simulación de los datos del ejemplo anterior, el diseño es balanceado (igual número de individuos en cada uno de los dos grupos y el mismo número de repeticiones en cada individuo de cada grupo en cada tratamiento), en mis datos reales es no - los dos grupos no tienen un número igual de individuos y el número de repeticiones varía. También, me imagino que el individuo debe ser establecido como un efecto aleatorio.
El trazado de la varianza de la muestra frente a la media de la muestra de la estimación de la probabilidad de éxito (denotado como p hat = K/N) de cada muestra de que la extrema probabilidades de éxito tienen menor varianza:
Este es eliminado cuando la estimación de probabilidades de éxito son transformados mediante la arcsen raíz cuadrada de la varianza de estabilización de transformación (denotado como arcsen(sqrt(p sombrero)):
El trazado de la varianza de la muestra de la transformación que calcula probabilidades de éxito frente a la media N muestra la esperada relación negativa:
El trazado de la varianza de la muestra de la transformación que calcula probabilidades de éxito para los dos grupos muestra que el grupo b tiene un poco más de varianzas, que es como yo simulada de los datos:
Finalmente, el trazado de la varianza de la muestra de la transformación que calcula probabilidades de éxito para los tres tratamientos no muestra ninguna diferencia entre los tratamientos, que es como yo simulada de los datos:
Hay alguna forma de un modelo lineal generalizado con el que puedo cuantificar el grupo y los efectos del tratamiento en la variación de las probabilidades de éxito?
Tal vez un heteroscedastic modelo lineal generalizado o alguna forma de loglineal de la varianza del modelo?
Algo en las líneas de un modelo en el que los modelos de la Varianza(y) = Zλ además de las E(y) = Xß, donde Z y X son los regresores de la media y la varianza, respectivamente, que en mi caso será idéntica y que incluyen el tratamiento (niveles de t.1, t.2, y t.3) y el grupo (niveles a y b), y, probablemente, N y R, y por lo tanto λ y β estimar sus efectos respectivos.
Como alternativa, podría ajustar un modelo a la muestra de varianzas a través de repeticiones de cada gen de cada grupo en cada tratamiento, utilizando un glm que los modelos que sólo tienen el valor esperado de la respuesta. La única pregunta aquí es cómo explicar el hecho de que diferentes genes tienen diferentes números de repeticiones. Creo que los pesos en un glm podría dar cuenta de que (muestra las desviaciones que se basan en más repeticiones deberían tener un mayor peso), pero exactamente lo que los pesos deben establecer?
Nota:
He intentado usar el dglm R paquete:
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
El efecto de grupo de acuerdo a dglm.el ajuste es bastante débil. Me pregunto si el modelo se establece el derecho o el poder de este modelo.
