如何绘制r中未知参数的幂和真实值的图形？

Question

假设IID随机示例X_I \ SIM BERNOULLI（\ pi） i = 1，\ dots，100和样本大小为100。我们想进行假设测试， H_0：\ pi \ ge 0.6，\，h_a：\ pi＆lt; 0.6 我想在我们的真实参数pi和电源之间获得关系。我有函数power。但是我只能输入每个TRUE pi = 0.50，0.51，0.52，0.53，...，0.59。如何绘制类似的数字？

我的代码如下。

n=100
pi0=0.60

##power function
power_fun=function(N,pi,alpha)
{
  pvalue=c()
  power=c()
  rej=c()
  
  for (i in 1:N) {
    set.seed(i)
    y=rbinom(n,size=1,prob=pi)
    pi_hat=mean(y)
    z=(pi_hat-pi0)/sqrt(1/n*pi_hat*(1-pi_hat))
    
    rej[i]=ifelse(z<qnorm(0.05,mean=0,sd=1),1,0)
    
    pvalue[i]=pnorm(q=z,mean=0,sd=1)
    
    c_value=qnorm(alpha,mean=0,sd=1)
    aug_term=(pi-pi0)/sqrt(1/n*pi_hat*(1-pi_hat))
    power[i]=pnorm(c_value-aug_term,mean=0,sd=1)
    
  }

  mean_pvalue=mean(pvalue)

  sd_pvalue=sd(pvalue)
  

  mean_power=ifelse(pi<pi0,mean(power),NA)

  sd_power=sd(power)
  
  rej_rate=mean(pvalue<alpha)
  
  type1_error=ifelse(pi>=pi0,rej_rate,NA)
  
  type2_error=ifelse(pi<pi0,1-mean_power,NA)
  
  type2_error_emp=ifelse(pi<pi0,1-rej_rate,NA)
  
  mc_se=sqrt(1/N*rej_rate*(1-rej_rate))
  
  df_out=data.frame("pvalue"=mean_pvalue,"pvalue2"=sd_pvalue,"power"=mean_power,"ESE_power"=sd_power, 
                    "rejection rate"=rej_rate,
                    "Type I error"=  type1_error,"type_II_error"=type2_error,"Type_II_error"=type2_error_emp,
                    "MC_SE"=mc_se)
  rownames(df_out)=paste0("N=",N,", pi=",pi, ", alpha=", alpha)
  
  df_out=round(df_out,3)
  
  return(df_out)
}

对于每个pi，我们可以获取电源。

power_fun(N=1000,pi=0.5,alpha=0.05)
#power=0.643
power_fun(N=1000,pi=0.51,alpha=0.05)
# power=0.566
power_fun(N=1000,pi=0.52,alpha=0.05)
power_fun(N=1000,pi=0.53,alpha=0.05)
power_fun(N=1000,pi=0.54,alpha=0.05)

但这太复杂了。是否有一种简单的方法来获取这些PI值的功能并绘制其pi和power的图表？

Answer 1

这更像是代码审查。

我认为您可以通过n，pi和alpha进行矢量化。
我认为您无法通过当前的实施。

遵循 @limey的建议：

#' Power function
#' 
#' 
#' @param N Total number of replications
#' @param pi 
#' @param alpha Significance level
#' @param pi0 
#' @param n Size of the replication
#'
#' @return
#' 
power_fun <- function(N, pi, alpha, pi0, n) {
  pvalue <- vector(mode = "numeric", length = N)
  power <- vector(mode = "numeric", length = N)
  rej <- vector(mode = "numeric", length = N)
  
  for (i in seq_len(N)) {
    # why?
    # set.seed(i)
    y = rbinom(n, size = 1, prob = pi)
    pi_hat = mean(y)
    z = (pi_hat - pi0) / sqrt(1 / n * pi_hat * (1 - pi_hat))
    
    rej[i] = ifelse(z < qnorm(0.05, mean = 0, sd = 1), 1, 0)
    
    pvalue[i] = pnorm(q = z, mean = 0, sd = 1)
    
    c_value = qnorm(alpha, mean = 0, sd = 1)
    aug_term = (pi - pi0) / sqrt(1 / n * pi_hat * (1 - pi_hat))
    power[i] = pnorm(c_value - aug_term, mean = 0, sd = 1)
    
  }
  
  mean_pvalue = mean(pvalue, na.rm = TRUE)
  sd_pvalue = sd(pvalue, na.rm = TRUE)
  # mean_power = ifelse(pi < pi0, mean(power), NA)
  mean_power = mean(power, na.rm = TRUE)
  sd_power = sd(power, na.rm = TRUE)
  rej_rate = mean(pvalue < alpha, na.rm = TRUE)
  
  type1_error = ifelse(pi >= pi0, rej_rate, NA)
  type2_error = ifelse(pi < pi0, 1 - mean_power, NA)
  type2_error_emp = ifelse(pi < pi0, 1 - rej_rate, NA)
  
  mc_se = sqrt(1 / N * rej_rate * (1 - rej_rate))
  
  df_out <- data.frame(
    N = N, 
    pi = pi, 
    alpha = alpha,
    pvalue = mean_pvalue,
    pvalue2 = sd_pvalue,
    power = mean_power,
    ESE_power = sd_power,
    rejection_rate = rej_rate,
    Type_I_error =  type1_error,
    type_II_error = type2_error,
    Type_II_error = type2_error_emp,
    MC_SE = mc_se
  )
  # moved this to above
  # rownames(df_out) = paste0("N=", N, ", pi=", pi, ", alpha=", alpha)
  
  # why?
  # df_out = round(df_out, 3)
  
  return(df_out)
}

请阅读评论和更改。

library(tidyverse)
n <- 100
pi0 <- 0.6

# testing the function
power_fun(N=1000,pi=0.5,alpha=0.05, pi0 = pi0, n = n)

#' Plot for all pi going from 0.5 to 1.
seq.default(0, 0.65, length.out = 200) %>%
# seq.default(0.5, 1., length.out = 2) %>% 
  map_df(function(pi) power_fun(N=1000,pi=pi,alpha=0.05, pi0 = pi0, n = n)) %>% 
  as_tibble() -> 
  power_df
# power_df %>% View()
power_df %>% {
  ggplot(., aes(pi, power)) + 
    geom_line() + 
    # geom_point() +
    geom_vline(xintercept = pi0, linetype = "dashed") +
    NULL
}