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 FiellerRatioCI
funkcji, 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