# # Priklad 96 # # Upozorneni: Nasledujici kod je napsan z duvodu pedagogicke nazornosti metody maxim. verohodnosti. # # Z hledika programovani v jazyku R by mnohem vyhodnejsi bylo pracovat maticove. # Data X1 <- c(1,2,2,2,3,3,3,4,4,4,4); X2 <- c(1,1,1,1,2,2,2,3,3,3,4); Y <- c(1,0,1,0,2,1,1,4,3,4,4) n <- length(X); # Log-verohodnosti loglik <- function(beta0, beta1, beta2) sum(Y*(beta0 + beta1*X1 + beta2*X2)) - sum(exp(beta0 + beta1*X1 + beta2*X2)) # # Skorova funkce # Prvni slozka U1 <- function(beta0, beta1, beta2) sum(Y) - sum(exp(beta0 + beta1*X1 + beta2*X2)) # Druha slozka U2 <- function(beta0, beta1, beta2) sum(Y*X1) - sum(exp(beta0 + beta1*X1 + beta2*X2)*X1) # Druha slozka U3 <- function(beta0, beta1, beta2) sum(Y*X2) - sum(exp(beta0 + beta1*X1 + beta2*X2)*X2) # MV odhad thetahat <- coef(glm(Y~X1+X2, family="poisson")) # Tato funkce umi najti MVO v tomto modelu (beta0hat <- thetahat[1]); # odhad parametru alfa (beta1hat <- thetahat[2]); # odhad parametru beta (beta2hat <- thetahat[3]); # odhad parametru beta # Kontrola, ze nase odhady opravdu resi skorove rovnice U1(beta0hat, beta1hat, beta2hat) # Mela by byt nula U2(beta0hat, beta1hat, beta2hat) # Mela by byt nula U3(beta0hat, beta1hat, beta2hat) # Mela by byt nula # # 1. Test pomerem verohodnosti # Odhad za H0 prum.Y <- mean(Y); beta0tilde <- log(prum.Y); # Testova statistika (LR <- 2*(loglik(beta0hat, beta1hat, beta2hat) - loglik(beta0tilde, 0, 0))) # p-hodnota (p.value.LR <- 1 - pchisq(LR, df=2)) # # 2. Rauv skorovy test # # Empiricka Fisherova informacni matice (FIM) J11 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)) J12 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)*X1) J13 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)*X2) J21 <- function(beta0, beta1, beta2) J12(beta0, beta1, beta2) J22 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)*X1^2) J23 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)*X1*X2) J31 <- function(beta0, beta1, beta2) J13(beta0, beta1, beta2) J32 <- function(beta0, beta1, beta2) J23(beta0, beta1, beta2) J33 <- function(beta0, beta1, beta2) 1/n*sum(exp(beta0 + beta1*X1 + beta2*X2)*X2^2) J <- function(beta0, beta1, beta2) matrix(c(J11(beta0, beta1, beta2), J21(beta0, beta1, beta2), J31(beta0, beta1, beta2), J12(beta0, beta1, beta2), J22(beta0, beta1, beta2), J23(beta0, beta1, beta2), J31(beta0, beta1, beta2), J32(beta0, beta1, beta2), J33(beta0, beta1, beta2)), 3, 3) # Inverze k empiricke FIM v bode thetatilde = (beta0tilde, 0) solve(J(beta0tilde, 0, 0)) # # My potrebujeme podmatici (2:3,2:3) vyse uvedene matice, tj. (Jinv2323 <- solve(J(beta0tilde, 0, 0))[2:3,2:3]); # Alternativne by sel pouzit i podmatice (2:3,2:3) inverze FIM v bode alfahat, betahat solve(J(beta0hat, beta1hat, beta2hat))[2:3,2:3] # Testova statistika Un <- c(U2(beta0tilde, 0, 0), U3(beta0tilde, 0, 0)); # potrebna cast skoru (Rn <- 1/n * t(Un)%*%Jinv2323%*%Un) # p-hodnota (p.value.Rn <- 1 - pchisq(Rn, df=1)) # # 3. Walduv test # Zde pouzijeme prvek (2,2) inverze FIM v bode alfahat, betahat (Jinv2323 <- solve(J(beta0hat, beta1hat, beta2hat))[2:3,2:3]); # Testova statistika (Wn <- n*t(c(beta1hat, beta2hat))%*%solve(Jinv2323)%*%c(beta1hat, beta2hat)); # p-hodnota (p.value.Wn <- 1 - pchisq(Wn, df=1)) # # Souhrn vysl.testovani <- cbind(c(LR, Rn, Wn), c(p.value.LR, p.value.Rn, p.value.Wn)); rownames(vysl.testovani) <- c("LR", "Rn", "Wn"); colnames(vysl.testovani) <- c("test. stat.", "p-hodnota"); print(vysl.testovani, digits=3) # # Intervalovy odhad o spolehlivosti 0.95 pro beta1 # (Jinv22 <- solve(J(beta0hat, beta1hat, beta2hat))[2,2]); beta1hat + c(-1,1) * qnorm(0.975) * sqrt(Jinv22)/sqrt(n)