R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("gap-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('gap') Loading required package: MASS Loading required package: genetics Loading required package: combinat Loading required package: gdata Loading required package: gtools Attaching package: 'genetics' The following object(s) are masked from package:MASS : genotype The following object(s) are masked from package:base : as.factor [1] "R/gap is loaded" Attaching package: 'gap' The following object(s) are masked from package:genetics : hap > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "bt" > > ### * bt > > flush(stderr()); flush(stdout()) > > ### Name: bt > ### Title: Bradley-Terry model for contingency table > ### Aliases: bt > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # Copeman JB, Cucca F, Hearne CM, Cornall RJ, Reed PW, > ##D # Ronningen KS, Undlien DE, Nistico L, Buzzetti R, Tosi R, et al. > ##D # (1995) Linkage disequilibrium mapping of a type 1 > ##D # diabetes susceptibility gene (IDDM7) to chromosome 2q31-q33. > ##D # Nat Genet 9: 80-5 > ##D > ##D x <- matrix(c(0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0, > ##D 0,0, 1, 3, 0,0, 0, 2, 3, 0, 0, 0, > ##D 2,3,26,35, 7,0, 2,10,11, 3, 4, 1, > ##D 2,3,22,26, 6,2, 4, 4,10, 2, 2, 0, > ##D 0,1, 7,10, 2,0, 0, 2, 2, 1, 1, 0, > ##D 0,0, 1, 4, 0,1, 0, 1, 0, 0, 0, 0, > ##D 0,2, 5, 4, 1,1, 0, 0, 0, 2, 0, 0, > ##D 0,0, 2, 6, 1,0, 2, 0, 2, 0, 0, 0, > ##D 0,3, 6,19, 6,0, 0, 2, 5, 3, 0, 0, > ##D 0,0, 3, 1, 1,0, 0, 0, 1, 0, 0, 0, > ##D 0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0, > ##D 0,0, 1, 0, 0,0, 0, 0, 0, 0, 0, 0),nrow=12) > ##D > ##D # Bradley-Terry model, only deviance is available in glm > ##D # (SAS gives score and Wald statistics as well) > ##D bt.ex<-bt(x) > ##D anova(bt.ex$bt.glm) > ##D summary(bt.ex$bt.glm) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "chow.test" > > ### * chow.test > > flush(stderr()); flush(stdout()) > > ### Name: chow.test > ### Title: Chow's test for heterogeneity in two regressions > ### Aliases: chow.test > ### Keywords: htest > > ### ** Examples > > ## Not run: > ##D dat1 <- matrix(c( > ##D 1.2, 1.9, 0.9, > ##D 1.6, 2.7, 1.3, > ##D 3.5, 3.7, 2.0, > ##D 4.0, 3.1, 1.8, > ##D 5.6, 3.5, 2.2, > ##D 5.7, 7.5, 3.5, > ##D 6.7, 1.2, 1.9, > ##D 7.5, 3.7, 2.7, > ##D 8.5, 0.6, 2.1, > ##D 9.7, 5.1, 3.6), byrow=TRUE, ncol=3) > ##D > ##D dat2 <- matrix(c( > ##D 1.4, 1.3, 0.5, > ##D 1.5, 2.3, 1.3, > ##D 3.1, 3.2, 2.5, > ##D 4.4, 3.6, 1.1, > ##D 5.1, 3.1, 2.8, > ##D 5.2, 7.3, 3.3, > ##D 6.5, 1.5, 1.3, > ##D 7.8, 3.2, 2.2, > ##D 8.1, 0.1, 2.8, > ##D 9.5, 5.6, 3.9), byrow=TRUE, ncol=3) > ##D > ##D y1<-dat1[,3] > ##D y2<-dat2[,3] > ##D x1<-dat1[,1:2] > ##D x2<-dat2[,1:2] > ##D chow.test.r<-chow.test(y1,x1,y2,x2) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "fbsize" > > ### * fbsize > > flush(stderr()); flush(stdout()) > > ### Name: fbsize > ### Title: Sample size for family-based linkage and association design > ### Aliases: fbsize > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D models <- matrix(c( > ##D 4.0, 0.01, > ##D 4.0, 0.10, > ##D 4.0, 0.50, > ##D 4.0, 0.80, > ##D 2.0, 0.01, > ##D 2.0, 0.10, > ##D 2.0, 0.50, > ##D 2.0, 0.80, > ##D 1.5, 0.01, > ##D 1.5, 0.10, > ##D 1.5, 0.50, > ##D 1.5, 0.80), ncol=2, byrow=TRUE) > ##D > ##D cat("\nThe family-based result: \n") > ##D cat("\ngamma p Y N_asp P_A Het N_tdt Het N_asp/tdt L_o L_s\n\n") > ##D for(i in 1:12) { > ##D g <- models[i,1] > ##D p <- models[i,2] > ##D fbsize(g,p) > ##D if(i%%4==0) cat("\n") > ##D } > ##D > ##D # APOE-4, Scott WK, Pericak-Vance, MA & Haines JL > ##D # Genetic analysis of complex diseases 1327 > ##D g <- 4.5 > ##D p <- 0.15 > ##D cat("\nAlzheimer's:\n\n") > ##D fbsize(g,p) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gc.em" > > ### * gc.em > > flush(stderr()); flush(stdout()) > > ### Name: gc.em > ### Title: Gene counting for haplotype analysis > ### Aliases: gc.em > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(hla) > ##D gc.em(hla[,3:8],locus.label=c("DQR","DQA","DQB")) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gcontrol" > > ### * gcontrol > > flush(stderr()); flush(stdout()) > > ### Name: gcontrol > ### Title: genomic control > ### Aliases: gcontrol > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D test<-c(1,2,3,4,5,6, 1,2,1,23,1,2, 100,1,2,12,1,1, > ##D 1,2,3,4,5,61, 1,2,11,23,1,2, 10,11,2,12,1,11) > ##D test<-matrix(test,nrow=6,byrow=T) > ##D gcontrol(test) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gcp" > > ### * gcp > > flush(stderr()); flush(stdout()) > > ### Name: gcp > ### Title: Permutation tests using GENECOUNTING > ### Aliases: gcp pm pmplus > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D > ##D data(fsnps) > ##D y<-fsnps$y > ##D cc<-1 > ##D g<-fsnps[,3:10] > ##D > ##D gcp(y,cc,g,miss.val="Z",n.sim=5) > ##D hap.score(y,g,method="hap",miss.val="Z") > ## End(Not run) > > > > cleanEx(); ..nameEx <- "genecounting" > > ### * genecounting > > flush(stderr()); flush(stdout()) > > ### Name: genecounting > ### Title: Gene counting for haplotype analysis > ### Aliases: genecounting > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # Now we use the HLA data for testing > ##D data(hla) > ##D hla.gc <- genecounting(hla[,3:8]) > ##D summary(hla.gc) > ##D hla.gc$l0 > ##D hla.gc$l1 > ##D > ##D # Now we use ALDH2 data > ##D data(aldh2) > ##D control <- gc.control(handle.miss=1) > ##D aldh2.gc <- genecounting(aldh2[,3:6],control=control) > ##D summary(aldh2.gc) > ##D aldh2.gc$l0 > ##D aldh2.gc$l1 > ##D > ##D # Now we perform gene counting with Chromosome X data > ##D # assuming allelic data have been extracted in columns 3-13 > ##D # and column 3 is sex > ##D dat <- mao[,3:13] > ##D loci <- c(12,9,6,5,3) > ##D contr <- gc.control(handle.miss=1) > ##D mao.gc <- genecounting(dat,xdata=TRUE,loci=loci,control=contr) > ##D mao.gc$npusr > ##D mao.gc$npdat > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gif" > > ### * gif > > flush(stderr()); flush(stdout()) > > ### Name: gif > ### Title: Kinship coefficient and genetic index of familiality > ### Aliases: gif > ### Keywords: datagen > > ### ** Examples > > ## Not run: > ##D test<-c( > ##D 5, 0, 0, > ##D 1, 0, 0, > ##D 9, 5, 1, > ##D 6, 0, 0, > ##D 10, 9, 6, > ##D 15, 9, 6, > ##D 21, 10, 15, > ##D 3, 0, 0, > ##D 18, 3, 15, > ##D 23, 21, 18, > ##D 2, 0, 0, > ##D 4, 0, 0, > ##D 7, 0, 0, > ##D 8, 4, 7, > ##D 11, 5, 8, > ##D 12, 9, 6, > ##D 13, 9, 6, > ##D 14, 5, 8, > ##D 16, 14, 6, > ##D 17, 10, 2, > ##D 19, 9, 11, > ##D 20, 10, 13, > ##D 22, 21, 20) > ##D test<-t(matrix(test,nrow=3)) > ##D gif(test,gifset=c(20,21,22)) > ##D > ##D # all individuals > ##D gif(test,gifset=1:23) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "hap" > > ### * hap > > flush(stderr()); flush(stdout()) > > ### Name: hap > ### Title: Haplotype reconstruction > ### Aliases: hap > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # 4 SNP example, to generate hap.out and assign.out alone > ##D data(fsnps) > ##D hap(id=fsnps[,1],data=fsnps[,3:10],nloci=4) > ##D dir() > ##D file.show("hap.out") > ##D file.show("assign.out") > ##D > ##D # to generate results of imputations > ##D control <- hap.control(ss=1,mi=5) > ##D hap(id=fsnps[,1],data=fsnps[,3:10],nloci=4,control=control) > ##D dir() > ## End(Not run) > > > > cleanEx(); ..nameEx <- "hap.em" > > ### * hap.em > > flush(stderr()); flush(stdout()) > > ### Name: hap.em > ### Title: Gene counting for haplotype analysis > ### Aliases: hap.em > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(hla) > ##D hap.em(id=1:length(hla[,1]),data=hla[,3:8],locus.label=c("DQR","DQA","DQB")) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "hap.score" > > ### * hap.score > > flush(stderr()); flush(stdout()) > > ### Name: hap.score > ### Title: Score Statistics for Association of Traits with Haplotypes > ### Aliases: hap.score > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D data(hla) > ##D y<-hla[,2] > ##D geno<-hla[,3:8] > ##D # complete data > ##D hap.score(y,geno,locus.label=c("DRB","DQA","DQB")) > ##D # incomplete genotype data > ##D hap.score(y,geno,locus.label=c("DRB","DQA","DQB"),handle.miss=1,mloci=1) > ##D unlink("assign.dat") > ##D > ##D ### note the differences in p values in the following runs > ##D data(aldh2) > ##D # to subset the data since hap doesn't handle one allele missing > ##D deleted<-c(40,239,256) > ##D aldh2[deleted,] > ##D aldh2<-aldh2[-deleted,] > ##D y<-aldh2[,2] > ##D geno<-aldh2[,3:18] > ##D # only one missing locus > ##D hap.score(y,geno,handle.miss=1,mloci=1,method="hap") > ##D # up to seven missing loci and with 10,000 permutations > ##D hap.score(y,geno,handle.miss=1,mloci=7,method="hap",n.sim=10000) > ##D > ##D # hap.score takes considerably longer time and does not handle missing data > ##D hap.score(y,geno,n.sim=10000) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "htr" > > ### * htr > > flush(stderr()); flush(stdout()) > > ### Name: htr > ### Title: Haplotype trend regression > ### Aliases: htr > ### Keywords: regression > > ### ** Examples > > ## Not run: > ##D # 26-10-03 > ##D # this is now part of demo > ##D test2<-read.table("test2.dat") > ##D y<-test2[,1] > ##D x<-test2[,-1] > ##D y<-as.matrix(y) > ##D x<-as.matrix(x) > ##D htr.test2<-htr(y,x) > ##D htr.test2 > ##D htr.test2<-htr(y,x,n.sim=10) > ##D htr.test2 > ##D > ##D # 13-11-2003 > ##D data(apoeapoc) > ##D apoeapoc.gc<-gc.em(apoeapoc[,5:8]) > ##D y<-apoeapoc$y > ##D for(i in 1:length(y)) if(y[i]==2) y[i]<-1 > ##D htr(y,apoeapoc.gc$htrtable) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "hwe" > > ### * hwe > > flush(stderr()); flush(stdout()) > > ### Name: hwe > ### Title: Hardy-Weinberg equlibrium test for multiallelic marker > ### Aliases: hwe > ### Keywords: htest > > ### ** Examples > > ## Not run: > ##D data(hla) > ##D hla.DQR <- hwe(hla[,3:4]) > ##D summary(hla.DQR) > ##D a <- c(3,2,2) > ##D a.out <- hwe(a,data.type="genotype") > ##D a.out > ##D a.out <- hwe(a,data.type="count") > ##D a.out > ## End(Not run) > > > > cleanEx(); ..nameEx <- "hwe.hardy" > > ### * hwe.hardy > > flush(stderr()); flush(stdout()) > > ### Name: hwe.hardy > ### Title: Hardy-Weinberg equilibrium test using MCMC > ### Aliases: hwe.hardy > ### Keywords: htest > > ### ** Examples > > ## Not run: > ##D # example 2 from hwe.doc: > ##D a<-c( > ##D 3, > ##D 4, 2, > ##D 2, 2, 2, > ##D 3, 3, 2, 1, > ##D 0, 1, 0, 0, 0, > ##D 0, 0, 0, 0, 0, 1, > ##D 0, 0, 1, 0, 0, 0, 0, > ##D 0, 0, 0, 2, 1, 0, 0, 0) > ##D ex2 <- hwe.hardy(a=a,alleles=8) > ##D > ##D # example using HLA > ##D data(hla) > ##D x <- hla[,3:4] > ##D y <- pgc(x,handle.miss=0,with.id=1) > ##D n.alleles <- max(x,na.rm=TRUE) > ##D z <- vector("numeric",n.alleles*(n.alleles+1)/2) > ##D z[y$idsave] <- y$wt > ##D hwe.hardy(a=z,alleles=n.alleles) > ##D > ##D # with use of class 'genotype' > ##D library(genetics) > ##D hlagen <- genotype(a1=x$DQR.a1, a2=x$DQR.a2, > ##D alleles=sort(unique(c(x$DQR.a1, x$DQR.a2)))) > ##D hwe.hardy(hlagen) > ##D > ##D # comparison with hwe > ##D hwe(z,data.type="count") > ##D > ##D # to create input file for HARDY > ##D print.tri<-function (xx,n) { > ##D cat(n,"\n") > ##D for(i in 1:n) { > ##D for(j in 1:i) { > ##D cat(xx[i,j]," ") > ##D } > ##D cat("\n") > ##D } > ##D cat("100 170 1000\n") > ##D } > ##D xx<-matrix(0,n.alleles,n.alleles) > ##D xxx<-lower.tri(xx,diag=TRUE) > ##D xx[xxx]<-z > ##D sink("z.dat") > ##D print.tri(xx,n.alleles) > ##D sink() > ##D # now call as: hwe z.dat z.out > ##D > ## End(Not run) > > > > cleanEx(); ..nameEx <- "kbyl" > > ### * kbyl > > flush(stderr()); flush(stdout()) > > ### Name: kbyl > ### Title: LD statistics for two multiallelic loci > ### Aliases: kbyl > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # example of two SNPs > ##D h <- c(0.442356,0.291532,0.245794,0.020319) > ##D n <- 481*2 > ##D kbyl(2,2,h,n) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "kin.morgan" > > ### * kin.morgan > > flush(stderr()); flush(stdout()) > > ### Name: kin.morgan > ### Title: kinship matrix for simple pedigree > ### Aliases: kin.morgan > ### Keywords: datagen > > ### ** Examples > > ## Not run: > ##D # Werner syndrome pedigree > ##D werner<-c( > ##D 1, 0, 0, 1, > ##D 2, 0, 0, 2, > ##D 3, 0, 0, 2, > ##D 4, 1, 2, 1, > ##D 5, 0, 0, 1, > ##D 6, 1, 2, 2, > ##D 7, 1, 2, 2, > ##D 8, 0, 0, 1, > ##D 9, 4, 3, 2, > ##D 10, 5, 6, 1, > ##D 11, 5, 6, 2, > ##D 12, 8, 7, 1, > ##D 13,10, 9, 2, > ##D 14,12, 11, 1, > ##D 15,14, 13, 1) > ##D werner<-t(matrix(werner,nrow=4)) > ##D kin.morgan(werner[,1:3]) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "makeped" > > ### * makeped > > flush(stderr()); flush(stdout()) > > ### Name: makeped > ### Title: A function to prepare pedigrees in post-MAKEPED format > ### Aliases: makeped > ### Keywords: datagen > > ### ** Examples > > ## Not run: > ##D library(gap) > ##D makeped("ped7.pre","ped7.ped",0,1,"ped7.lop") > ## End(Not run) > > > > cleanEx(); ..nameEx <- "mia" > > ### * mia > > flush(stderr()); flush(stdout()) > > ### Name: mia > ### Title: multiple imputation analysis for hap > ### Aliases: mia > ### Keywords: utilities > > ### ** Examples > > ## Not run: > ##D # 4 SNP example, to generate hap.out and assign.out alone > ##D data(fsnps) > ##D hap(id=fsnps[,1],gdata=fsnps[,3:10],nloci=4) > ##D > ##D # to generate results of imputations > ##D hap(id=fsnps[,1],gdata=fsnps[,3:10],nloci=4,ss=1,mi=5) > ##D > ##D # to extract information from the second run above > ##D mia(so=1,ns=1,mi=5) > ##D file.show("mia.out") > ##D > ##D ## commands to check out where the output files are as follows: > ##D ## Windows > ##D # system("command.com") > ##D ## Unix > ##D # system("csh") > ## End(Not run) > > > > cleanEx(); ..nameEx <- "mtdt" > > ### * mtdt > > flush(stderr()); flush(stdout()) > > ### Name: mtdt > ### Title: Transmission/disequilibrium test of a multiallelic marker > ### Aliases: mtdt > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # Copeman et al (1995) Nat Genet 9: 80-5 > ##D > ##D x <- matrix(c(0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0, > ##D 0,0, 1, 3, 0,0, 0, 2, 3, 0, 0, 0, > ##D 2,3,26,35, 7,0, 2,10,11, 3, 4, 1, > ##D 2,3,22,26, 6,2, 4, 4,10, 2, 2, 0, > ##D 0,1, 7,10, 2,0, 0, 2, 2, 1, 1, 0, > ##D 0,0, 1, 4, 0,1, 0, 1, 0, 0, 0, 0, > ##D 0,2, 5, 4, 1,1, 0, 0, 0, 2, 0, 0, > ##D 0,0, 2, 6, 1,0, 2, 0, 2, 0, 0, 0, > ##D 0,3, 6,19, 6,0, 0, 2, 5, 3, 0, 0, > ##D 0,0, 3, 1, 1,0, 0, 0, 1, 0, 0, 0, > ##D 0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0, > ##D 0,0, 1, 0, 0,0, 0, 0, 0, 0, 0, 0),nrow=12) > ##D > ##D # See note to bt for the score test obtained by SAS > ##D > ##D mtdt(x) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "muvar" > > ### * muvar > > flush(stderr()); flush(stdout()) > > ### Name: muvar > ### Title: Means and variances under 1- and 2- locus (biallelic) QTL model > ### Aliases: muvar > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # the default 1-locus model > ##D muvar(n.loci=1,y1=c(0,1,1),p1=0.5) > ##D > ##D # the default 2-locus model > ##D muvar(n.loci=2,y12=c(1,1,1,1,1,0,0,0,0),p1=0.99,p2=0.9) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "pbsize" > > ### * pbsize > > flush(stderr()); flush(stdout()) > > ### Name: pbsize > ### Title: Power for population-based association design > ### Aliases: pbsize > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D models <- matrix(c( > ##D 4.0, 0.01, > ##D 4.0, 0.10, > ##D 4.0, 0.50, > ##D 4.0, 0.80, > ##D 2.0, 0.01, > ##D 2.0, 0.10, > ##D 2.0, 0.50, > ##D 2.0, 0.80, > ##D 1.5, 0.01, > ##D 1.5, 0.10, > ##D 1.5, 0.50, > ##D 1.5, 0.80), ncol=2, byrow=TRUE) > ##D > ##D g <- 4.5 > ##D p <- 0.15 > ##D cat("\nAlzheimer's:\n\n") > ##D > ##D zalpha <- 5.45 # 5.4513104 > ##D z1beta <- -0.84 > ##D > ##D q <- 1-p > ##D pi <- 0.065 # 0.07 generates 163, equivalent to ASP > ##D k <- pi*(g*p+q)^2 > ##D s <- (1-pi*g^2)*p^2+(1-pi*g)*2*p*q+(1-pi)*q^2 > ##D # LGL formula > ##D lambda <- pi*(g^2*p+q-(g*p+q)^2)/(1-pi*(g*p+q)^2) > ##D # my own > ##D lambda <- pi*p*q*(g-1)^2/(1-pi*(g*p+q)^2) > ##D # not sure about +/-! > ##D n <- (z1beta+zalpha)^2/lambda > ##D > ##D # may be used to correct for population prevalence > ##D cat("\nThe population-based result: Kp=",k, "Kq=",s, "n=",ceiling(n),"\n") > ##D > ##D # population-based sample size > ##D strlen <- function(x) length(unlist(strsplit(as.character(x),split=""))) > ##D kp <- c(0.01,0.05,0.10) > ##D cat("\nRandom ascertainment with disease prevalence\n") > ##D cat("\n 1% 5% 10%\n\n") > ##D for(i in 1:12) { > ##D g <- models[i,1] > ##D p <- models[i,2] > ##D q <- 1-p > ##D for(j in 1:3) { > ##D n <- pbsize(g,p,kp[j]) > ##D cat(rep("",12-strlen(ceiling(n))),format(ceiling(n))) > ##D } > ##D cat("\n") > ##D if(i%%4==0) cat("\n") > ##D } > ##D cat("This is only an approximation, a more accurate result\n") > ##D cat("can be obtained by Fisher's exact test\n") > ## End(Not run) > > > > cleanEx(); ..nameEx <- "pedtodot" > > ### * pedtodot > > flush(stderr()); flush(stdout()) > > ### Name: pedtodot > ### Title: Converting pedigree(s) to dot file(s) > ### Aliases: pedtodot > ### Keywords: dplot > > ### ** Examples > ## Not run: > ##D # > ##D # Three examples of pedigree-drawing > ##D # assuming pre-MakePed LINKAGE file in which IDs are characters > ##D pre<-read.table("pheno.pre",as.is=T)[,1:6] > ##D pedtodot(pre) > ##D dir() > ##D # for post-MakePed LINKAGE file in which IDs are integers > ##D ped <-read.table("pheno.ped")[,1:10] > ##D pedtodot(ped,makeped=T) > ##D dir() > ##D # for a single file with a list of pedigrees ordered data > ##D sink("gaw14.dot") > ##D pedtodot(ped,sink=F) > ##D sink() > ##D file.show("gaw14.dot") > ##D # more details > ##D pedtodot(ped,sink=F,page="B5",url="http://www.ucl.ac.uk/~rmjdjhz") > ## End(Not run) > > > > cleanEx(); ..nameEx <- "pfc" > > ### * pfc > > flush(stderr()); flush(stdout()) > > ### Name: pfc > ### Title: Probability of familial clustering of disease > ### Aliases: pfc > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # IPF among 203 siblings of 100 COPD patients from Liang KY, SL Zeger, B Qaquish (1992) > ##D # Multivariate regression analyses for categorical data (with discussion). J Roy Stat Soc > ##D # B 54:3-40 > ##D > ##D # the degrees of freedom is 15 > ##D famtest<-c( > ##D 1, 0, 36, > ##D 1, 1, 12, > ##D 2, 0, 15, > ##D 2, 1, 7, > ##D 2, 2, 1, > ##D 3, 0, 5, > ##D 3, 1, 7, > ##D 3, 2, 3, > ##D 3, 3, 2, > ##D 4, 0, 3, > ##D 4, 1, 3, > ##D 4, 2, 1, > ##D 6, 0, 1, > ##D 6, 2, 1, > ##D 6, 3, 1, > ##D 6, 4, 1, > ##D 6, 6, 1) > ##D test<-t(matrix(famtest,nrow=3)) > ##D famp<-pfc(test) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "pfc.sim" > > ### * pfc.sim > > flush(stderr()); flush(stdout()) > > ### Name: pfc.sim > ### Title: Probability of familial clustering of disease > ### Aliases: pfc.sim > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D > ##D # Li, Fraumeni, Mulvihill, et al., 1988 > ##D # family_size #_of_affected frequency > ##D > ##D famtest<-c( > ##D 1, 0, 2, > ##D 1, 1, 0, > ##D 2, 0, 1, > ##D 2, 1, 4, > ##D 2, 2, 3, > ##D 3, 0, 0, > ##D 3, 1, 2, > ##D 3, 2, 1, > ##D 3, 3, 1, > ##D 4, 0, 0, > ##D 4, 1, 2, > ##D 5, 0, 0, > ##D 5, 1, 1, > ##D 6, 0, 0, > ##D 6, 1, 1, > ##D 7, 0, 0, > ##D 7, 1, 1, > ##D 8, 0, 0, > ##D 8, 1, 1, > ##D 8, 2, 1, > ##D 8, 3, 1, > ##D 9, 3, 1) > ##D > ##D test<-matrix(famtest,byrow=T,ncol=3) > ##D > ##D famp<-pfc.sim(test) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "pgc" > > ### * pgc > > flush(stderr()); flush(stdout()) > > ### Name: pgc > ### Title: Preparing weight for GENECOUNTING > ### Aliases: pgc > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D > ##D data(hla) > ##D x <- hla[,3:8] > ##D > ##D # do not handle missing data > ##D y<-pgc(x,handle.miss=0,with.id=1) > ##D hla.gc<-genecounting(y$cdata,y$wt,handle.miss=0) > ##D > ##D # handle missing but with multilocus genotype identifier > ##D pgc(x,handle.miss=1,with.id=1) > ##D > ##D # handle missing data with no identifier > ##D pgc(x,handle.miss=1,with.id=0) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "plot.hap.score" > > ### * plot.hap.score > > flush(stderr()); flush(stdout()) > > ### Name: plot.hap.score > ### Title: Plot Haplotype Frequencies versus Haplotype Score Statistics > ### Aliases: plot.hap.score > ### Keywords: misc > > ### ** Examples > > ## Not run: > ##D save <- hap.score(y, geno, trait.type = "gaussian") > ##D > ##D # Example illustrating generic plot function: > ##D plot(save) > ##D > ##D # Example illustrating specific method plot function: > ##D plot.hap.score(save) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "print.hap.score" > > ### * print.hap.score > > flush(stderr()); flush(stdout()) > > ### Name: print.hap.score > ### Title: Print a hap.score object > ### Aliases: print.hap.score > ### Keywords: print > > ### ** Examples > > ## Not run: > ##D save <- hap.score(y, geno, trait.type = "gaussian") > ##D > ##D # Example illustrating generic print function: > ##D print(save) > ##D > ##D # Example illustrating specific method print function: > ##D print.hap.score(save) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "s2k" > > ### * s2k > > flush(stderr()); flush(stdout()) > > ### Name: s2k > ### Title: Statistics for 2 by K table > ### Aliases: s2k > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D # an example from Mike Neale > ##D # termed 'ugly' contingency table by Patrick Sullivan > ##D y1 <- c(2,15,16,35,132,30,25,7,12,24,10,10,0) > ##D y2 <- c(0, 6,31,49,120,27,15,8,14,25, 3, 9,3) > ##D > ##D result <- s2k(y1,y2) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "tbyt" > > ### * tbyt > > flush(stderr()); flush(stdout()) > > ### Name: tbyt > ### Title: LD statistics for two SNPs > ### Aliases: tbyt > ### Keywords: models > > ### ** Examples > > ## Not run: > ##D h <- c(0.442356,0.291532,0.245794,0.020319) > ##D n <- 481*2 > ##D tbyt(h,n) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "whscore" > > ### * whscore > > flush(stderr()); flush(stdout()) > > ### Name: whscore > ### Title: Whittemore-Halpern scores for allele-sharing > ### Aliases: whscore > ### Keywords: utilities > > ### ** Examples > > ## Not run: > ##D c<-matrix(c(1,1,1,2,2,2),ncol=2) > ##D whscore(c,type=1) > ##D whscore(c,type=2) > ## End(Not run) > > > > ### *