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.