R help - Jun 2011 - Hardy Weinberg Simulation

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]]

> 
> 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.
What do you expect? Depending on the genetic model, you may not see HWE in the
"cases".
> datamat[h,] <- t(rmultinom(1, size=c(10, 40, 50), prob=c(0.33, 0.33,
0.34)))
"size" should be the total size ie 100.

-- 
| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v