Chapter 4 Hierarchical Simulations

This simulation study aims to deal with the second goal specified in the introduction. After an effect has been discovered, how can data from replication studies be combined with the original? In this scenario, we must account for the heterogeneity between sites; neglecting the uncertainty that comes from replication studies would lead to erroneously confident estimates of significance and effect size.

4.1 Data Generation Procedure

The data are generated from a hierarchical (i.e. mixed effect) logistic model as discribed in the models section: if truly associated, \(\mu\) and \(\sigma^2\) have (fixed) nonzero values;\(\beta_j \sim \textsf{N}(\mu, \sigma^2)\). Otherwise, \(\mu=\beta_j = 0, \forall j\).

To try to keep this simulation as close to the real data as possible, a preliminary logistic regression with random slope and random p53 coefficient by site was run. This led to the values of \(\mu =0.203 , \sigma^2 = 0.003\). The value of \(\mu\) remained fixed through all the simulations, but different values of \(\sigma\) were used to test the sensitivity of the models: \(\sigma\), \(\sigma/2\), \(2\sigma\), and \(4\sigma\). The number of sites was set to 7, since results using 30 sites were almost identical. Each site had 1000 observations. A total of 100 simulated datasets was created each time.

100 datasets were also simulated under the null hypothesis. They were fit with the models described previously.

4.1.1 Finding “Discovery” Sites

In this simulation study, a logistic regression with fixed effects for sites was conducted to find the site with the smallest p-value less than \(\alpha\). If no sites matched this description, the data was resampled until at least one site was viable. The maximum likelihood estimate of this effect and its variance were added as data for the conditional likelihood model, and the p-value was added to the bayes factor approximation model. The observations for this site were then taken out of the data.

4.2 Results

4.2.1 RMSE of \(\mu\)

As expected, the models performed more poorly as \(\sigma\) increased. Out of the three proposed models (compared with the original), the fully bayesian and BF approximation models performed best when there was no true effect (since they were the only ones that that had this option). However, there were some simulated datasets where the bayes factor model estimate was actually nonzero and quite large, which suggests that it is not nearly as reliable as the bayesian model.

At small variances (s/2, s, 2s), the original and bayesian models outperform the others. This is not surprising since the other models only have access to \(\frac{6}{7}\) of the data. The Bayesian model actually has a slightly higher lower RMSE than the other models when there is a true association.

Table 4.1: RMSE of \(\mu\)

	mu=0	s/2	s	2s	4s
Bayesian	0.000	0.025	0.037	0.053	0.099
Cond Likelihood	0.018	0.027	0.034	0.055	0.098
BF Approx	0.012	0.024	0.031	0.053	0.097
Original	0.024	0.025	0.031	0.052	0.099

4.2.2 Coverage of \(\mu\)

The conditional likelihood and BF approximation models are the most conservative, with the intervals covering 0 more times than the others for large values of \(\sigma\). All models have very high coverage in general.

The Bayesian model had the shortest intervals, and the original model had the largest. Thus, even though the coverage and RMSE are around the same, the new models are preferable to the original. This does not apply to the simulation with \(\sigma=4s\), because \(4s>\mu\), which leads to more negative site effects. Thus, it makes sense for models to have wider credible intervals for these simulations.

4.2.3 Probability of Association \(\xi\)

While one would expect the probability of being associated (\(\xi\)) to also increase with \(\sigma\), this was not true for either the fully bayesian model nor the bayes factor approximation one, both of which had consistent posterior estimates of \(\xi\). Similarly, the proportion of nonzero \(\mu\) samples from the posterior (this is the same as the proportion of times the latent variable \(i = 1\)), was almost 1 for the truly associated cases, and close to zero for true null. One thing to consider is that under the null hypothesis, the variance across sites would actually be zero, which is why the models identified the association so decisively.

For the simulations that had no true effect, the Bayes Factor approximation model has much larger median \(\xi\) and greater proportion of sampled \(H_1\) because the mass of the distribution of \(\xi\) is shifted towards 1 by the Bayes Factor transformation. Thus, even though the proportion of \(H_1\) is quite low (and the median of \(\mu\) is 0), the mean of \(\xi\) is greater than \(0.5\).

Table 4.2: Average of Posterior Median \(\xi\)

	mu=0	s/2	s	2s	4s
Bayesian	0.210	0.836	0.830	0.832	0.837
BF Approx	0.847	0.991	0.992	0.994	0.998

Table 4.3: Average Proportion of \(H_1\)

	mu=0	s/2	s	2s	4s
Bayesian	0.132	0.998	0.989	0.989	0.996
BF Approx	0.501	1.000	1.000	1.000	1.000

4.2.4 Sensitivity of Conditional Likelihood Method to Changes in \(\alpha\)

The conditional likelihood method with random effects is robust to changes in the level \(\alpha\). To test this, we consider 5 different levels: \(0.05, 0.01, 0.005, 0.001, 10^{-7}\). 100 datasets were sampled, for which at least one location was significant at the smallest \(\alpha\) level. The conditional likelihood model with random effects and without was then fitted for each level.

Table 4.4: RMSE of \(\mu\)

	0	s	4s
0.05	0.019	0.039	0.101
0.01	0.020	0.040	0.090
0.005	0.019	0.040	0.099
0.001	0.019	0.038	0.098
1e-07	NA	0.039	0.097

This model shows little difference across levels of \(\alpha\).