Jadi tujuan akhir saya adalah memiliki plot dengan beberapa interval kepercayaan 95% yang diplot secara vertikal dalam 2 grup seperti contoh ini:

Text

Saya telah menemukan kode ini: https://rpubs.com/nati147/703463

Tapi bagaimana saya bisa menambahkan perbandingan groupwise dalam plot?

Tulis fungsi 'CI_95' yang mengambil input vektor nilai sampel, dan menghasilkan interval kepercayaan 95% untuk sampel ini. Anda dapat menggunakan fungsi 'margin_error_95'.

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

Tulis fungsi yang disebut 'margin_error_95' yang mengambil input vektor nilai sampel, dan menampilkan margin kesalahan untuk interval kepercayaan 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
}

Jalankan plot:

plot_CI_95(1)

Text

0
H. berg 12 Mei 2021, 16:25

1 menjawab

Jawaban Terbaik

Berikut adalah solusi, meskipun mungkin memerlukan beberapa penyempurnaan lebih lanjut berdasarkan preferensi Anda. Saya mempertahankan struktur umum fungsi plot_CI_95 Anda, tetapi menambahkan loop pada grup yang berbeda. Ini berarti bahwa variabel mu dan sig sekarang harus memiliki beberapa nilai, satu untuk setiap grup jika Anda ingin menunjukkan perbedaan berdasarkan grup. Ada juga beberapa warna dan parameter grafis lainnya. Hasilnya ditunjukkan di bawah ini.

Untuk menghindari interval untuk dua grup yang tumpang tindih, ada beberapa parameter yang harus diubah. 1) angka height di png (meningkatkan nilai membuat lebih banyak ruang antar grup, 2) parameter offset (dapat meningkat hingga sekitar 0,3), atau 3) lwd di segments fungsi (nilai yang lebih kecil berarti garis yang lebih tipis). Menggunakan png atau fungsi serupa untuk menyimpan gambar secara langsung akan memungkinkan penyempurnaan tampilan yang diinginkan.

Example Plot with Group Differences

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()
1
Calvin 17 Mei 2021, 13:59