垂直 95% 置信区间图 2 组比较

Question

所以我的最终目标是在 2 个组中垂直绘制多个 95% 置信区间的图，如下例所示：

我找到了这个代码：https://rpubs.com/nati147/703463

但是如何在图中添加分组比较？

编写一个函数“CI_95”，输入一个样本值向量， 并输出该样本的 95% 置信区间。您可以使用 “margin_error_95”函数。

CI_95 <- function(sample_vals, sig){
  error <- margin_error_95(sample_vals, sig)
  CI <- mean(sample_vals) + c(-1, 1)*error
}

编写一个名为“margin_error_95”的函数，它接受输入一个向量 样本值，并输出 95% 置信度的误差范围 间隔。

margin_error_95 <- function(sample_vals, sig){
  n <- length(sample_vals)
  mar_err <- 1.96*(sig/sqrt(n))
}

plot_CI_95 <- function(seed){
  B <- 100
  n <- 30
  mu <- 5
  sig <- 1.2
  
  set.seed(seed)
  # extract upper bound of CI's
  
  x_1 <- replicate(B,
                   {samp <- rnorm(n, mean = mu, sd = sig )
                   max(CI_95(samp, sig))
                   }
  )
  
  #extract lower bound of CI's
  
  set.seed(seed)
  
  x_0 <- replicate(B,
                   {samp <- rnorm(n, mean = mu, sd = sig )
                   min(CI_95(samp, sig))
                   }
  )
  
  set.seed(seed)
  
  means <- replicate(B, mean(rnorm(n, mean = mu, sd = sig)))
  
  plot(means, 1:B, pch = 20,
       xlim = c(mu - sig, mu + sig),
       ylim = c(0,B+1),
       xlab = "sample means",
       ylab = "index of the CI",
       main = paste(B, "Confidence intervals")
  )
  
  for (i in 1:B){
    if(between(mu, x_0[i], x_1[i])){
      segments(x_0[i], i, x_1[i], i, lwd = 2) #plot CI's that contain the mean in black
    } else {
      segments(x_0[i], i, x_1[i], i, col = "red", lwd = 2) #plot CI's that don't contain the mean in red
    }
  }
  
  abline(v=mu, col = "blue") #plot a vertical line at the population mean
}

运行剧情：

plot_CI_95(1)

Answer 1

这里有一个解决方案，尽管它可能需要根据您的偏好进行一些进一步的改进。我保留了 plot_CI_95 函数的一般结构，但在不同的组上添加了一个循环。这意味着mu 和sig 变量现在必须有多个值，如果要显示组间差异，则每个组一个值。还有一些颜色和其他图形参数。结果如下所示。

为避免两组的间隔重叠，需要调整一些参数。 1）height in png（增加值使组之间的空间更大，2）offset 参数（可以增加到大约 0.3），或 3）lwd 在segments 函数中（较小的值意味着细线）。使用png 或类似函数直接保存图形将允许微调所需的外观。

library(dplyr)

CI_95 <- function(sample_vals, sig){
  error <- margin_error_95(sample_vals, sig)
  CI <- mean(sample_vals) + c(-1, 1)*error
}


margin_error_95 <- function(sample_vals, sig){
  n <- length(sample_vals)
  mar_err <- 1.96*(sig/sqrt(n))
}

png("group_plot.png",height=7,width=3,units = 'in',res=1000)

plot_CI_95 <- function(seed){
  B <- 100
  n <- 30
  # mean and std dev as a vector of values
  # need to have same length
  # assume one value per group
  mu <- c(5,4.5) # group1, group2
  sig <- c(1.2,1) # group1, group2
  offset<- 0.25 # controls point and line offset from nominal value
  # colors for 2 groups
  colsuse  <- c('steelblue','gold')
  
  # loop over groups
  # mu and sig are now indexed by this loop
  for(j in 1:length(mu)){
    # use seed+j to make different random sample for each group
    
    # extract upper bound of CI's
    set.seed(seed+j)
    x_1 <- replicate(B,
                     {samp <- rnorm(n, mean = mu[j], sd = sig[j])
                     max(CI_95(samp, sig[j]))
                     }
    )
    
    #extract lower bound of CI's
    set.seed(seed+j)
    x_0 <- replicate(B,
                     {samp <- rnorm(n, mean = mu[j], sd = sig[j])
                     min(CI_95(samp, sig[j]))
                     }
    )
    
    set.seed(seed+j)
    
    means <- replicate(B, mean(rnorm(n, mean = mu[j], sd = sig[j])))
    
    # for first group, establish the plot
    # for second group, add values to the plot
    # if groups are very different this might need to be modified with the xlim
    if(j == 1){
      plot(means, (1:B)+offset*ifelse(j==1,1,-1), pch = 20,
           xlim = c(mu[j] - sig[j], mu[j] + sig[j]),
           ylim = c(0,B+1),
           xlab = "sample means",
           ylab = "index of the CI",
           main = paste(B, "Confidence intervals"),
           col=colsuse[j]
      )
    }else{
      points(means, (1:B)+offset*ifelse(j==1,1,-1), pch = 20,
           col=colsuse[j])
    }
 
    for (i in 1:B){
      if(between(mu[j], x_0[i], x_1[i])){
        segments(x_0[i], i+offset*ifelse(j==1,1,-1), x_1[i], i+offset*ifelse(j==1,1,-1), col=colsuse[j], lwd = 1) #plot CI's that contain the mean in black
      } else {
        segments(x_0[i], i+offset*ifelse(j==1,1,-1), x_1[i], i+offset*ifelse(j==1,1,-1), col = "red", lwd = 1) #plot CI's that don't contain the mean in red
      }
    }
    
    abline(v=mu[j], col = colsuse[j]) #plot a vertical line at the population mean
  }
  
  par(xpd=F)
  # legend for the groups
  legend("topright",legend = c('Male','Female'),lty=1,col=colsuse,cex=0.5)
}

plot_CI_95(1)

dev.off()