Jak obliczyć przedział ufności stosunku dwóch normalnych średnich

26

Chcę ustalić limity dla przedziału ufności dla stosunku dwóch średnich. Załóżmy, że X 1N ( θ 1 , σ 2 ) i X 2N ( θ 2 , σ 2 ) są niezależne, a średni stosunek Γ = θ 1 / θ 2 . Próbowałem rozwiązać: Pr ( - z ( α / 2100(1-α)%
X1N.(θ1,σ2))X2)N.(θ2),σ2))Γ=θ1/θ2)ale równania tego nie można rozwiązać w wielu przypadkach (bez pierwiastków). czy robię coś źle? Czy istnieje lepsze podejście? Dzięki

Par(-z(α/2)))X1-ΓX2)/σ1+γ2)z(α/2)))=1-α
francogrex
źródło
1
Problem polega na tym, że stosunek dwóch liczb z dwóch rozkładów normalnych jest zgodny z rozkładem Cauchy'ego, a zatem wariancja jest niezdefiniowana.
6
@mbq - rozkład Cauchy'ego nie przedstawia problemów dla przedziałów ufności, ponieważ CDF jest odwrotną funkcją styczną. Wariancja nie musi być zdefiniowana, aby elementy CI mogły działać. A stosunek dwóch normalnych RV ze średnią zerową wynosi Cauchy, ale niekoniecznie dwóch normalnych RV ze średnią niezerową.
prawdopodobieństwo prawdopodobieństwa
@probabilityislogic Na pewno muszę przestać myśleć w niedzielne poranki.

Odpowiedzi:

31

Metoda Fiellera robi to, co chcesz - oblicza przedział ufności dla ilorazu dwóch średnich, przy założeniu, że oba są próbkowane z rozkładów Gaussa.

  • Oryginalny cytat to: Fieller EC: Biologiczna standaryzacja insuliny. Dostarcz do JR Statist Soc 1940, 7: 1-64.

  • Artykuł w Wikipedii dobrze podsumowuje.

  • Stworzyłem kalkulator online, który wykonuje obliczenia.

  • Oto strona podsumowująca matematykę z pierwszego wydania mojej intuicyjnej biostatystyki

Harvey Motulsky
źródło
To bardzo dobre referencje, podoba mi się również, że zrobiłeś dla niego kalkulator (+1). Jak można się jednak spodziewać, w kalkulatorze wyraźnie stwierdza się, że gdy przedział ufności mianownika zawiera zero, nie można obliczyć CI ilorazu. Myślę, że to samo dzieje się, gdy próbuję rozwiązać równanie kwadratowe. załóżmy, że wariancja wynosi 1, mu1 = 0, a mu2 = 1, N = 10000. To nierozwiązywalne.
francogrex
2
dzięki za kalkulatora online Harvey, jestem typowym biologiem z niewystarczającym zapleczem statystycznym, a kalkulator był dokładnie tym, czego potrzebowałem.
Timtico,
Niesamowity kalkulator - dokładnie to, czego szukałem. Dzięki
Alexander
@ harvey-motulsky link do dodatku już nie działa. Zastanawiałem się, czy materiał z tego dodatku jest zawarty w trzecim wydaniu Intuitive Biostatistics?
Gabriel Południowy
@GabrielSouthern Dziękujemy za wskazanie zgnilizny linków. Naprawiony.
Harvey Motulsky
1

Również jeśli chcesz obliczyć przedział ufności Fiellera, który nie używa mratios(zwykle dlatego, że nie chcesz prostego dopasowania lm, ale na przykład dopasowania glmer lub glmer.nb), możesz użyć następującej FiellerRatioCIfunkcji, modelując dane wyjściowe modelu, aname nazwa parametru licznika, b name nazwa parametru denomiatora. Można także użyć bezpośrednio funkcji FiellerRatioCI_basic, która daje a, b oraz macierz kowariancji między a i b.

Zauważ, że alfa ma tutaj wartość 0,05 i jest „zakodowane” na 1,96 w kodzie. Możesz je zastąpić dowolnym poziomem Studenta.

FiellerRatioCI <- function (x, ...) { # generic Biomass Equilibrium Level
    UseMethod("FiellerRatioCI", x)
}
FiellerRatioCI_basic <- function(a,b,V,alpha=0.05){
    theta <- a/b
    v11 <- V[1,1]
    v12 <- V[1,2]
    v22 <- V[2,2]

    z <- qnorm(1-alpha/2)
    g <- z*v22/b^2
    C <- sqrt(v11 - 2*theta*v12 + theta^2 * v22 - g*(v11-v12^2/v22))
    minS <- (1/(1-g))*(theta- g*v12/v22 - z/b * C)
    maxS <- (1/(1-g))*(theta- g*v12/v22 + z/b * C)
    return(c(ratio=theta,min=minS,max=maxS))
}
FiellerRatioCI.glmerMod <- function(model,aname,bname){
    V <- vcov(model)
    a<-as.numeric(unique(coef(model)$culture[aname]))
    b<-as.numeric(unique(coef(model)$culture[bname]))
    return(FiellerRatioCI_basic(a,b,V[c(aname,bname),c(aname,bname)]))
}
FiellerRatioCI.glm <- function(model,aname,bname){
    V <- vcov(model)
    a <- coef(model)[aname]
    b <- coef(model)[bname]
    return(FiellerRatioCI_basic(a,b,V[c(aname,bname),c(aname,bname)]))
}

Przykład (na podstawie standardowego podstawowego przykładu glm):

 counts <- c(18,17,15,20,10,20,25,13,12)
 outcome <- gl(3,1,9)
 treatment <- gl(3,3)
 glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

 FiellerRatioCI(glm.D93,"outcome2","outcome3")
ratio.outcome2            min            max 
      1.550427      -2.226870      17.880574
cmbarbu
źródło