How to simulate phenotype from two random multivariant normal distributed variable and calculate statistical power!

Question

Greeting Homo Sapiens,

I am a Biology student and just stepped into the area of Statistical Genomics. I have gone through several posts from the forum and I must say the quality response from the community helped me to solve my 20-30% problem. Still I am finding some gap in my research area for which I am trying to write this question. Lets hope I can explain my problem and make you guys understand easily 0_0

So basically the aim of my study is to

1. Perform Simulation through Mixed linear model
2. Calculate Power of Statistical Models
3. Find Type 1 error in my model

Now my model is a bit tricky (at least for me since I am from Bioscience department). I am trying to fit a model with Genotypic Variance, two random variables and these variables should be multivariant normally distributed. Next I am going to select some fixed positions from Genotype variance variable and multiply with two random multivariant normal distributed variable to simulate Phenotype (This will help me with finding Phenotypic Variance explained by my model). Finally I want to calculate the Power and Type 1 error of my model. The steps are mentioned below:

1. Create Genotype Variance from Genotype .raw file
2. Select any three positions from the variance variable
3. Create First random effect by using mvnorm function from MASS package
4. Create Second random effect by using mvnorm function
5. Combine first three random effect columns with variance variable
6. Create Simulation effect from two random variables (First three columns)
7. Combine the simulated effect with variance variables i.e. effect_indiv_geno<-t(GEN_QTN)%*%effect_1_simu %*% effect_2_simu
8. Create residual from rnorm
9. Calculate Phenotype from effect_indiv_geno and residual

I need help from here onwards as I am stuck. I am not sure how to plan the Power calculation process and then follow type 1 error for the model. Looking forward to hearing from you guys. Thanks and Stay Safe.

Answer 1

It's difficult (impossible?) to calculate the a-priori power of a mixed model with formulas. The best approach is to simulate not only your data but also your sampling from the data.

Simulate a very large data set via your steps so far, with your intended assignments of Phenotypes to the combinations of random-variable values. That provides the "alternative hypothesis" for your study of power/Type-II error. Then make a copy of that data set, but shuffle the Phenotype assignments randomly. That provides a null condition: no (intended) associations of Phenotype with the random variables. The null-condition data allow you to evaluate Type I-error.

Then sample (with replacement) the number of cases that you would intend to take in your actual study, from each simulated data set. Do this repeatedly, say 999 times. Run your model separately on each sample from each of those data sets. For each of those 2 data sets, keep all the coefficient estimates based on each sample. When you're done, for each of those 2 data sets, examine the empirical distributions of the 999 values of each of the coefficient estimates returned by your model.

Type I error, by definition, is the fraction of cases showing a "significant" result when there really is none--what you find on modeling the samples of null-condition data. If you want the limits for p < 0.05 "significance" based on the null, for each coefficient those would be the 25th and 975th of your 999 values (extreme 2.5% on each end of the distribution).

You also can use a nominal 5% probability cutoff with the tools provided by a package like lmer test on your 999 samples from the null, and see whether you actually find only 5% false positives. That way, you can see how well the nominal significance cutoff you specify for your model matches up with how it works in practice, under the null condition.

Type-II error (1 - power) is the fraction of times you miss a true positive result, based on an assumed cutoff for Type-I error: what you find on modeling samples from the "alternative hypothesis" data set. As you designed that data set to have a true relationship between the random-variable values and Phenotype, what fraction of the 999 models erroneously show no relationship based on your cutoff for Type-I error?

With this type of sampling, you can directly examine the inherent tradeoffs between actual Type I and Type II errors as you adjust the nominal significance cutoff. By changing the size of the samples taken from the large data set, you can estimate the sample size that you would need to obtain desired levels of Type I and Type II errors.

That's the general outline. This page has some hints for simulating mixed models. Also, be aware that the calculation or interpretation of p-values from mixed models isn't straightforward; see this page for an introduction and links.