gap (version 1.1-20)

hwe.hardy: Hardy-Weinberg equilibrium test using MCMC

Description

Hardy-Weinberg equilibrium test by MCMC

Usage

hwe.hardy(a, alleles = 3, seed = 3000, sample = c(1000, 1000, 5000))

Arguments

a

an array containing the genotype counts, as integer.

alleles

number of allele at the locus, greater than or equal to 3, as integer

seed

pseudo-random number seed, as integer.

sample

optional, parameters for MCMC containing number of chunks, size of a chunk and burn-in steps, as integer.

Value

The returned value is a list containing:

method

Hardy-Weinberg equilibrium test using MCMC

data.name

name of used data if x is given

p.value

Monte Carlo p value

p.value.se

standard error of Monte Carlo p value

switches

percentage of switches (partial, full and altogether)

References

Guo, S.-W. and E. A. Thompson (1992) Performing the exact test of Hardy-Weinberg proportion for multiple alleles. Biometrics. 48:361--372.

See Also

hwe, HWE.test, genotype

Examples

# NOT RUN {
  
# }
# NOT RUN {
    # example 2 from hwe.doc:
    a<-c(
    3,
    4, 2,
    2, 2, 2,
    3, 3, 2, 1,
    0, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 1,
    0, 0, 1, 0, 0, 0, 0,
    0, 0, 0, 2, 1, 0, 0, 0)
    ex2 <- hwe.hardy(a=a,alleles=8)

    # example using HLA
    data(hla)
    x <- hla[,3:4]
    y <- pgc(x,handle.miss=0,with.id=1)
    n.alleles <- max(x,na.rm=TRUE)
    z <- vector("numeric",n.alleles*(n.alleles+1)/2)
    z[y$idsave] <- y$wt
    hwe.hardy(a=z,alleles=n.alleles)

    # with use of class 'genotype'
    # this is to be fixed
    library(genetics)
    hlagen <- genotype(a1=x$DQR.a1, a2=x$DQR.a2, 
                       alleles=sort(unique(c(x$DQR.a1, x$DQR.a2))))
    hwe.hardy(hlagen)

    # comparison with hwe
    hwe(z,data.type="count")

    # to create input file for HARDY
    print.tri<-function (xx,n) {
        cat(n,"\n")
        for(i in 1:n) {
            for(j in 1:i) {
                cat(xx[i,j]," ")
            }
        cat("\n")
        }
        cat("100 170 1000\n")
    }
    xx<-matrix(0,n.alleles,n.alleles)
    xxx<-lower.tri(xx,diag=TRUE)
    xx[xxx]<-z
    sink("z.dat")
    print.tri(xx,n.alleles)
    sink()
    # now call as: hwe z.dat z.out
  
# }