# Score and Wald tests (also LRT and Asymp. LRT) library(latex2exp) rm(list=ls()) ########################################################################### ########################################################################### ##### Pareto distribution ########################################################################### ########################################################################### beta <- seq(2.5,6.5,length = 41) beta0 <- 4 alpha <- 0.05 n <- 100 S <- 5000 powSc <- powALR <- powZ1 <- powZ2 <- powLR <- numeric(length(beta)) # find the critical values of the LR test: rej when T < c1 or T > c2, where T = beta0/beta.hat m <- 1000 d <- numeric(m) u <- c(1:m) / (m + 1) * alpha for( k in 1:m){ c1 <- qgamma(u[k],n,scale = 1/n) c2 <- qgamma(1 - alpha + u[k],n,scale = 1/n) d[k] <- abs(c1 * exp(-c1) - c2 * exp(-c2)) } c1 <- qgamma(u[which.min(d)],n,scale = 1/n) c2 <- qgamma(1 - alpha + u[which.min(d)],n,scale = 1/n) # compute the exact power curve for the LR test. powLR.theo <- pgamma(beta/beta0 * c1,n,scale = 1/n) + 1 - pgamma(beta/beta0 * c2,n,scale = 1/n) for(j in 1:length(beta)){ Sc <- Z1 <- Z2 <- ALR <- LR <- numeric(S) for( s in 1:S){ U <- runif(n) X <- exp( - (1/beta[j]) * log(1 - U) ) beta.hat <- 1/mean(log(X)) # Score test Sc[s] <- n * beta0^2 *(1/beta.hat - 1/beta0)^2 # Wald test 1 Z1[s] <- sqrt(n) * ( beta.hat - beta0) / beta.hat # Wald test 2 Z2[s] <- sqrt(n) * ( beta.hat - beta0) / beta0 # Asymptotic LRT ALR[s] <- - 2 * n * (log(beta0 / beta.hat) - (beta0/beta.hat - 1)) # LR test (complicated to find the critical values) LR[s] <- beta0/beta.hat } rejSc <- Sc > qchisq(1-alpha,1) rejZ1 <- abs(Z1) > qnorm(1-alpha/2) rejZ2 <- abs(Z2) > qnorm(1-alpha/2) rejALR <- ALR > qchisq(1-alpha,1) rejLR <- (LR < c1) | (LR > c2) powSc[j] <- mean(rejSc) powZ1[j] <- mean(rejZ1) powZ2[j] <- mean(rejZ2) powALR[j] <- mean(rejALR) powLR[j] <- mean(rejLR) } plot(NA, ylim = c(0,1), xlim = range(beta), xlab = "beta", ylab = paste("power ( n = ",n,", beta0 = ",beta0," )", sep = "")) lines(powSc ~ beta, col = "chartreuse", lwd = 5) lines(powALR ~ beta, col = "blue", lwd = 1) lines(powZ1 ~ beta) lines(powZ2 ~ beta, lty = 2) # lines(powLR ~ beta, lty = 1, col = "blue") lines(powLR.theo ~ beta, lty = 4, col = "dark red") abline( h = alpha, lty = 3) abline( v = beta0, lty = 3) x.pos <- grconvertX(.7, from = "nfc", to = "user") y.pos <- grconvertY(.7, from = "nfc", to = "user") legend(x = x.pos, y = y.pos, legend = c("Score", "Asymp LR", "Wald 1", "Wald 2","LR"), col = c("chartreuse", "blue", "black", "black","dark red"), lwd = c(5,1,1,1,1), lty = c(1,1,1,2,4), bty = "n") #### Now compare the values of the statistics: j <- 21 Sc <- Z1 <- Z2 <- ALR <- LR <- numeric(S) for( s in 1:S){ U <- runif(n) X <- exp( - (1/beta[j]) * log(1 - U) ) beta.hat <- 1/mean(log(X)) # Score test Sc[s] <- n * beta0^2 *(1/beta.hat - 1/beta0)^2 # Wald test 1 Z1[s] <- sqrt(n) * ( beta.hat - beta0) / beta.hat # Wald test 2 Z2[s] <- sqrt(n) * ( beta.hat - beta0) / beta0 # Asymptotic LRT ALR[s] <- - 2 * n * (log(beta0 / beta.hat) - (beta0/beta.hat - 1)) # LR test (complicated to find the critical values) LR[s] <- beta0/beta.hat } cbind(Z1,Z2,ALR,Sc,LR)[1:15,] ########################################################################### ########################################################################### ##### Test hypotheses about the Weibull shape parameter with score test ########################################################################### ########################################################################### rm(list=ls()) n <- 30 a0 <- 1.5 a <- seq(.5,3,length = 31) b <- 2 alpha <- 0.05 S <- 5000 powSc <- numeric(length(a)) for(j in 1:length(a)){ rejS <- numeric(S) for( s in 1:S){ U <- runif(n) X <- b*(log(1/(1-U)))^(1/a[j]) b0.hat <- (mean(X^a0))^(1/a0) # Score test (advantage is you do not have to compute the unrestricted MLE): S0 <- matrix(c(n/a0 + sum(log(X/b0.hat)) - sum( (X/b0.hat)^a0 * log(X/b0.hat)), -n*a0/b0.hat + a0 / b0.hat * sum( (X/b0.hat)^a0 )),nrow = 2) I0 <- n * matrix(c(1.824/a0^2,-0.423/b0.hat, -0.423/b0.hat,a0^2/b0.hat^2 ), nrow = 2, ncol = 2, byrow = TRUE) Sstat <- t(S0) %*% solve(I0) %*% S0 rejS[s] <- Sstat > qchisq(1-alpha,1) } powSc[j] <- mean(rejS) } # plot power curves plot(NA, ylim = c(0,1), xlim = range(a), xlab = "a", ylab = paste("power of Score test ( n = ",n,", a0 = ",a0,", b = ",b," )",sep= "")) lines(powSc ~ a) abline(h = alpha, lty = 3) abline(v = a0, lty = 3) ########################################################################### ########################################################################### ##### Test hypotheses about the bernoulli proportion with various tests ########################################################################### ########################################################################### rm(list = ls()) n <- 20 p0 <- 1/3 p <- seq(0.01,0.99, length = 99) alpha <- 0.05 S <- 10000 powSc <- powZ1 <- powZ2 <- powW <- numeric(length(p)) for( j in 1:length(p)){ # generate S data sets Y <- rbinom(S, n, p[j]) p.hat <- Y / n # Score statistic Sstat <- n * ( p.hat * (1 - p0) - (1 - p.hat)*p0)^2/(p0*(1-p0)) rejSc <- Sstat > qchisq(1 - alpha, 1) powSc[j] <- mean(rejSc) # Wald test 1 Z1 <- sqrt(n) * (p.hat - p0) / sqrt(p.hat *(1-p.hat)) rejZ1 <- abs(Z1) > qnorm(1 - alpha/2) powZ1[j] <- mean(rejZ1) # Wald test 2 Z2 <- sqrt(n) * (p.hat - p0) / sqrt(p0 *(1-p0)) rejZ2 <- abs(Z2) > qnorm(1 - alpha/2) powZ2[j] <- mean(rejZ2) # asymptotic LRT W <- - 2 * n * ( p.hat * log( p0/p.hat) + (1 - p.hat) * log((1 - p0)/(1 - p.hat))) rejW <- W > qchisq(1 - alpha,1) powW[j] <- mean(rejW) } # plot power curves plot(NA, ylim = c(0,1), xlim = c(0,1), xlab = "p", ylab = paste("power ( n = ",n,", p0 = ",round(p0,3)," )",sep = "")) lines(powSc ~ p, col = "chartreuse", lwd = 5) lines(powW ~ p, col = "blue", lwd = 3) lines(powZ1 ~ p) lines(powZ2 ~ p, lty = 2) abline(h = alpha, lty = 3) abline(v = p0, lty = 3) legend(x = .7, y = .8, legend = c("Score", "Asymp LR", "Wald 1", "Wald 2"), col = c("chartreuse", "blue", "black", "black"), lwd = c(5,3,1,1), lty = c(1,1,1,2), bty = "n") # Now make a table showing some of the test statistic values, to show that they are pretty similar: j <- 25 # generate S data sets Y <- rbinom(S, n, p[j]) p.hat <- Y / n # Score statistic Sstat <- n * ( p.hat * (1 - p0) - (1 - p.hat)*p0)^2/(p0*(1-p0)) rejSc <- Sstat > qchisq(1 - alpha, 1) powSc[j] <- mean(rejSc) # Wald Z1 <- sqrt(n) * (p.hat - p0) / sqrt(p.hat *(1-p.hat)) rejZ1 <- abs(Z1) > qnorm(1 - alpha/2) powZ1[j] <- mean(rejZ1) Z2 <- sqrt(n) * (p.hat - p0) / sqrt(p0 *(1-p0)) rejZ2 <- abs(Z2) > qnorm(1 - alpha/2) powZ2[j] <- mean(rejZ2) # asymptotic LRT W <- - 2 * n * ( p.hat * log( p0/p.hat) + (1 - p.hat) * log((1 - p0)/(1 - p.hat))) rejW <- W > qchisq(1 - alpha,1) powW[j] <- mean(rejW) cbind(Z1,Z2,W,Sstat)[1:15,] #### Now do a size comparison of the tests for different choice of p0 rm(list = ls()) n <- 20 p0.seq <- seq(0.01,0.99, length = 99) alpha <- 0.05 S <- 10000 sizeSc <- sizeZ1 <- sizeZ2 <- sizeW <- numeric(length(p0.seq)) for( j in 1:length(p0.seq)){ p0 <- p0.seq[j] # generate S data sets Y <- rbinom(S, n, p0) p.hat <- Y / n # Score statistic Sstat <- n * ( p.hat * (1 - p0) - (1 - p.hat)*p0)^2/(p0*(1-p0)) rejSc <- Sstat > qchisq(1 - alpha, 1) sizeSc[j] <- mean(rejSc) # Wald test 1 Z1 <- sqrt(n) * (p.hat - p0) / sqrt(p.hat *(1-p.hat)) rejZ1 <- abs(Z1) > qnorm(1 - alpha/2) sizeZ1[j] <- mean(rejZ1) # Wald test 2 Z2 <- sqrt(n) * (p.hat - p0) / sqrt(p0 *(1-p0)) rejZ2 <- abs(Z2) > qnorm(1 - alpha/2) sizeZ2[j] <- mean(rejZ2) # asymptotic LRT W <- - 2 * n * ( p.hat * log( p0/p.hat) + (1 - p.hat) * log((1 - p0)/(1 - p.hat))) rejW <- W > qchisq(1 - alpha,1) sizeW[j] <- mean(rejW) } # plot power curves plot(NA, ylim = c(0,.9), xlim = c(0,1), xlab = "p0", ylab = paste("size ( n = ",n," )",sep = "")) lines(sizeSc ~ p0.seq, col = "chartreuse", lwd = 5) lines(sizeW ~ p0.seq, col = "blue", lwd = 3) lines(sizeZ1 ~ p0.seq) lines(sizeZ2 ~ p0.seq, lty = 2) abline(h = alpha, lty = 3) legend(x = .7, y = .7, legend = c("Score", "Asymp LR", "Wald 1", "Wald 2"), col = c("chartreuse", "blue", "black", "black"), lwd = c(5,3,1,1), lty = c(1,1,1,2), bty = "n")