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ę:
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)):
Wykreślenie wariancji próbki przekształconych oszacowanych prawdopodobieństw sukcesu w porównaniu ze średnią N pokazuje oczekiwaną ujemną zależność:
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:
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:
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.