---
title: "yeast power study"
author: "H Qin"
date: "8/31/2021"
output:
pdf_document: default
html_document: default
---
```{r simulate genotypes}
rm(list=ls())
N = 150
nuc_means = rpois(10, lambda=10)
mit_means = rpois(15, lambda=10)
summary(nuc_means)
summary(mit_means)
b0= 0
b1= 1 # mito influence on phenotype
b2= 1 # nuclear influence on phenotype
b3 =0.2 # mit X nuc interaction influence on phenotype, p << 0.001
b3 =0.1 # mit X nuc interaction influence on phenotype, p=0.049
#b3 = 0.05 # p = 0.3
```
```{r simulate-phenotype}
debug = 0
phenotype_mit_nuc = function(b0, b1, b2, b3, mit_single_mean, nuc_single_mean, debug){
y = b0 + b1*mit_single_mean + b2*nuc_single_mean + b3*mit_single_mean * nuc_single_mean
if (debug > 0) {
print( paste("pmn:: mit_single_mean =", mit_single_mean, "nuc_single_mean", nuc_single_mean) )
}
return (y)
}
nuc_genotypes = sample(1:10, N, replace=TRUE)
mit_genotypes = sample(1:15, N, replace=TRUE)
y = 1:N
for ( i in 1:N ){
#print(paste("i:", i, "mit_genotypes[i]",mit_genotypes[i] ))
y[i] = phenotype_mit_nuc(b0, b1, b2, b3, mit_means[mit_genotypes[i]], nuc_means[nuc_genotypes[i]], debug=0) + rnorm(1)
}
tb = data.frame( cbind( y, mit_genotypes, nuc_genotypes))
tb$mit_genotypes = factor( tb$mit_genotypes)
tb$nuc_genotypes = factor( tb$nuc_genotypes)
summary(tb)
```
```{r}
library(nlme);
m1a = glm(y ~ mit_genotypes , tb, family='gaussian');
m2 = glm(y ~ mit_genotypes + nuc_genotypes , tb, family='gaussian');
m3 = glm( y ~ mit_genotypes + nuc_genotypes + mit_genotypes:nuc_genotypes, data=tb)
```
```{r}
#summary(m1a)
```
```{r}
anova( m1a, m2, test='F')
```
```{r}
summary(m2)
summary(m3)
anova(m2, m3, test='F')
```
