Simulations/simulate_meta_data.function.R at master · agaiasaracini/Simulations · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134


#To simulate one data set of a meta analysis with and without ORB
simulate_meta_data <- function(n_studies, mu, tau_squared, n_treatment=50, n_control=50, gamma){


  #Number of studies per meta-analysis
  n_studies <- n_studies

  #True treatment effect and true heterogeneity parameter
  mu <- mu
  tau_squared <- tau_squared

  #Empty data frame
  meta_data <- data.frame(study = 1:n_studies,
                          n_control = numeric(n_studies),
                          n_treatment = numeric(n_studies),
                          y = numeric(n_studies),
                          s = numeric(n_studies),
                          p_value = numeric(n_studies))


  n_treatment <- rep(n_treatment, n_studies)
  n_control <- rep(n_control, n_studies)


  # Generate true treatment effect mean and standard errors for all studies at once
  theta_i <- rnorm(n_studies, mean = mu, sd = sqrt(tau_squared))
  sigma_i <- sqrt(rchisq(n_studies, df = 2*n_control - 2) / ((n_control - 1) * n_control))

  # Generate observed treatment effect for all studies at once
  y_i <- rnorm(n_studies, mean = theta_i, sd = sigma_i)

  # Calculate one-sided p-values for all studies at once
  p_i <- pnorm(y_i / sigma_i, lower.tail = FALSE)
  p_i <- round(p_i, 10)


  # Fill in data frame
  meta_data$n_control <- n_control
  meta_data$n_treatment <- n_treatment
  meta_data$y <- y_i
  meta_data$s <- sigma_i
  meta_data$p_value <- p_i


  #Missing data mechanism
  meta_data_miss <- meta_data


  #probability of being missing, ie 1-prob of selection
  replace_prob <- exp(-4 *as.numeric(meta_data$p_value)^gamma)
  high_index <- runif(n_studies) >= replace_prob
  meta_data_miss$y[high_index] <- "high"
  meta_data_miss$s[high_index] <- "high"

  #Probability of being unreported and HR
  #replace_prob <- exp(-4 *as.numeric(meta_data$p_value)^gamma)

  #for (j in 1:n_studies){

   # if (runif(1) >= replace_prob[j]) {
    #  meta_data_miss$y[j] <- "high"
     # meta_data_miss$s[j] <- "high"
  #  } else {
   #   meta_data_miss$y[j] <- meta_data_miss$y[j]
    #  meta_data_miss$s[j] <- meta_data_miss$s[j]
    #}

  #}


  if (length(which(meta_data_miss$y == "high")) > 0){
    non_high_indices <- which(meta_data_miss$y != "high")
    n_high_rows <- length(which(meta_data_miss$y == "high"))
    n_values_to_replace <- min(sample(0:n_high_rows, 1), length(non_high_indices))

    if (length(non_high_indices) > 1){
      replace_indices <- sample(non_high_indices, n_values_to_replace)
    } else {
      replace_indices <- non_high_indices
    }

    #replace_indices <- sample(non_high_indices, n_values_to_replace)
    meta_data_miss$y[replace_indices] <- "low"
    meta_data_miss$s[replace_indices] <- "low"
  } else {
    meta_data_miss$y <- meta_data_miss$y
    meta_data_miss$s <- meta_data_miss$s
  }


  #among the rest, choose some at random some to be low risk
  #if (length(which(meta_data_miss$y == "high")) > 0) { # if there is at least one HR
   # non_high_indices <- which(meta_data_miss$y != "high") # indices reported
  #  n_high_rows <- length(which(meta_data_miss$y == "high")) # number of HR
   # n_values_to_replace <- min(sample(1:n_high_rows, 1), length(non_high_indices)) # number of values to replace so that never more low than high
  #  replace_indices <- ifelse(length(non_high_indices) == 1,
   #                           non_high_indices, #if there is only one non-high study, pick that one to be low
    #                          sample(non_high_indices, n_values_to_replace)) #if there are more than one, pick
  #  meta_data_miss$y[replace_indices] <- "low"
  #  meta_data_miss$s[replace_indices] <- "low"

  #} else {

  #  meta_data_miss$y <- meta_data_miss$y
   # meta_data_miss$s <- meta_data_miss$s

  #}

  #make sure we have at least 3 studies that are reported otheriwse calculationg of ci unstable
  #pls also v unrealistic that there would be a meta analysis conducted on 2 studies!
  non_high_low_rows <- which(!(meta_data_miss$y %in% c("high", "low")))
  n_non_high_low_rows <- length(non_high_low_rows)

  if (n_non_high_low_rows < 3) {
    high_low_rows <- which(meta_data_miss$y %in% c("high", "low"))

    if (n_non_high_low_rows == 0) {
      random_rows <- sample(high_low_rows, 3)
    } else {
      random_rows <- sample(high_low_rows, 3 - n_non_high_low_rows)
    }

    meta_data_miss$y[random_rows] <- meta_data$y[random_rows]
    meta_data_miss$s[random_rows] <- meta_data$s[random_rows]
  }

  meta_data_miss$s_true <- meta_data$s

  return(meta_data_miss)

}