Hello,
I am trying to simulate 10 relicates of 100-tables. Each table is a 2 x 3
and 80% pf the tables are true nulls and 20% are non-nulls. The nulls follow
the Hardy Weinberg distribution (ratio) 1:2:1.
I have the code below but the p-values are not what I am expecting. I want
to use the Cochran Armitage trend test to get the p-values.
num.reps=10
num.vars=1000
pi0 = 80
num.subjects = 100
#Create list control list.control and list.treatment
#Minor alleles will be from (0.10, 0.30)
#num.subjects is the number of subjects in each
datamat = matrix(0,num.vars,3)
list.treat <- vector("list", 1) # create list
for (i in 1: num.reps)
{
list.treat[[i]] <- matrix(0, nrow = num.vars, ncol=3)
}
for ( k in 1: num.reps)
{
minor.allele.prop = runif(num.vars, 0.10, 0.30)
###################################
##### Generating under H1 #####
###################################
## Minor Allele Frequency ##
p.control <- minor.allele.prop # p = P(a) #
q.control <- 1-p.control # q = P(A) = 1-P(a) #
## Genotype by Mendelian
rule: P(G) ##
p.aa.control <- p.control^2 # P(G=aa) #
p.Aa.control <- 2*p.control*q.control # P(G=Aa) #
p.AA.control <- q.control^2 # P(G=AA) #
###################################
##### Generating Control Data #####
###################################
## Minor Allele Frequency ##
for (h in 1:800)
{
datamat[h,] <- t(rmultinom(1, size=c(10, 40, 50), prob=c(0.33, 0.33, 0.34)))
}
for (h in 801:1000)
{
datamat[h,] <- t(rmultinom(1, size=c(10,40,50), prob=c(p.aa.control[h],
p.Aa.control[h], p.AA.control[h])))
}
list.treat[[k]] <- datamat
}
--
Thanks,
Jim.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
