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.

> ### * <HEADER>
> ###
> 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("ldDesign-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('ldDesign')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "SS.oneway.bf"
> 
> ### * SS.oneway.bf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SS.oneway.bf
> ### Title: Bayes factors for one-way analysis of variance models.
> ### Aliases: SS.oneway.bf
> ### Keywords: htest models
> 
> ### ** Examples
> 
> # Bayes factors corresponding to P-values 0.05,0.01,0.001,0.0001 for n=200
> SS.oneway.bf(group.sizes=c(50,100,50),Fstat=qf(0.95,2,197))
[1] 0.4186708
> SS.oneway.bf(group.sizes=c(50,100,50),Fstat=qf(0.99,2,197))
[1] 2.145294
> SS.oneway.bf(group.sizes=c(50,100,50),Fstat=qf(0.999,2,197))
[1] 22.21853
> SS.oneway.bf(group.sizes=c(50,100,50),Fstat=qf(0.9999,2,197))
[1] 230.1143
> 
> 
> 
> cleanEx(); ..nameEx <- "ld.design"
> 
> ### * ld.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ld.design
> ### Title: Functions for design of experiments to detect linkage
> ###   disequilibrium
> ### Aliases: ld.design ld.power
> ### Keywords: design
> 
> ### ** Examples
> 
> ld.power(n=seq(100,1000,by=100),p=0.5,q=0.5,D=0.1,h2=0.1,phi=0,Bf=20)
         n      power
 [1,]  100 0.02202122
 [2,]  200 0.04173622
 [3,]  300 0.07258736
 [4,]  400 0.11431801
 [5,]  500 0.16610790
 [6,]  600 0.22646948
 [7,]  700 0.29304168
 [8,]  800 0.36345961
 [9,]  900 0.43512021
[10,] 1000 0.50576803
attr(,"parms")
             p              q              D             h2            phi 
           0.5            0.5            0.1            0.1            0.0 
            Bf missclass.rate 
          20.0            0.0 
> ld.design(p=0.5,q=0.5,D=0.1,h2=0.1,phi=0,Bf=20,power=0.9,print.it=TRUE,nmin=600,nmax=4000)
              n     power
 [1,]  600.0000 0.2264695
 [2,]  623.6856 0.2417236
 [3,]  648.3062 0.2580209
 [4,]  673.8987 0.2751798
 [5,]  700.5016 0.2933890
 [6,]  728.1546 0.3126167
 [7,]  756.8992 0.3327295
 [8,]  786.7785 0.3540140
 [9,]  817.8374 0.3762007
[10,]  850.1223 0.3992875
[11,]  883.6818 0.4234433
[12,]  918.5660 0.4483356
[13,]  954.8273 0.4740780
[14,]  992.5200 0.5005715
[15,] 1031.7008 0.5275421
[16,] 1072.4282 0.5551778
[17,] 1114.7633 0.5830947
[18,] 1158.7697 0.6111580
[19,] 1204.5133 0.6394231
[20,] 1252.0627 0.6674401
[21,] 1301.4891 0.6952016
[22,] 1352.8667 0.7225033
[23,] 1406.2724 0.7490147
[24,] 1461.7864 0.7748022
[25,] 1519.4919 0.7994795
[26,] 1579.4754 0.8229126
[27,] 1641.8267 0.8450837
[28,] 1706.6395 0.8656927
[29,] 1774.0107 0.8847592
[30,] 1844.0416 0.9021798
[31,] 1916.8370 0.9178525
[32,] 1992.5060 0.9318985
[33,] 2071.1622 0.9442272
[34,] 2152.9233 0.9549255
[35,] 2237.9121 0.9641075
[36,] 2326.2559 0.9718240
[37,] 2418.0872 0.9782410
[38,] 2513.5436 0.9834741
[39,] 2612.7682 0.9876581
[40,] 2715.9099 0.9909647
[41,] 2823.1231 0.9935080
[42,] 2934.5687 0.9954307
[43,] 3050.4137 0.9968549
[44,] 3170.8318 0.9978820
[45,] 3296.0036 0.9986090
[46,] 3426.1166 0.9991089
[47,] 3561.3659 0.9994437
[48,] 3701.9544 0.9996628
[49,] 3848.0927 0.9998013
[50,] 4000.0000 0.9998865
attr(,"parms")
             p              q              D             h2            phi 
           0.5            0.5            0.1            0.1            0.0 
            Bf missclass.rate 
          20.0            0.0 
[1] 1834.544
> ld.design(p=0.5,q=0.5,D=0.1,h2=0.1,phi=0,Bf=20,power=0.9,print.it=FALSE,nmin=1700,nmax=1900)
[1] 1834.724
> 
> 
> 
> cleanEx(); ..nameEx <- "ld.sim"
> 
> ### * ld.sim
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ld.sim
> ### Title: Functions to simulate populations with a bi-allelic marker and
> ###   QTL in linkage disequilibrium and test for association.
> ### Aliases: ld.sim ld.sim1
> ### Keywords: design
> 
> ### ** Examples
> 
> # Power from stochastic simulation for Luo's population 12
> data(luo.ld.populations)
> luo.pop12.sim <- ld.sim(nsim=3000,
+                         n=luo.ld.populations[12,"n"],
+                         p=luo.ld.populations[12,"p"],
+                         q=luo.ld.populations[12,"q"],
+                         D=luo.ld.populations[12,"D"],
+                         h2=luo.ld.populations[12,"h2"],
+                         phi=luo.ld.populations[12,"phi"],
+                         Vp=100)
...............................
> # power
> table(luo.pop12.sim[,4] < 0.05)[2]/sum(table(luo.pop12.sim[,4] < 0.05))
     TRUE 
0.5393548 
> # Cf power from deterministic calculation
> luo.ld.power(n=luo.ld.populations[12,"n"],
+              p=luo.ld.populations[12,"p"],
+              q=luo.ld.populations[12,"q"],
+              D=luo.ld.populations[12,"D"],
+              h2=luo.ld.populations[12,"h2"],
+              phi=luo.ld.populations[12,"phi"],
+              Vp=100,
+              alpha=0.05)
ems.beta = 370.91 , ems.omega = 97.25 , power = 0.5454 
> 
> 
> 
> cleanEx(); ..nameEx <- "luo.ld.populations"
> 
> ### * luo.ld.populations
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: luo.ld.populations
> ### Title: Luo's Linkage disequilibrium example populations
> ### Aliases: luo.ld.populations
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(luo.ld.populations)
> luo.ld.populations
      pop.   n   p   q   D  h2 phi
 [1,]    1 100 0.5 0.5 0.1 0.1 0.0
 [2,]    2 200 0.5 0.5 0.1 0.1 0.0
 [3,]    3 200 0.5 0.5 0.2 0.1 0.0
 [4,]    4 200 0.5 0.5 0.1 0.2 0.0
 [5,]    5 200 0.5 0.5 0.1 0.1 0.5
 [6,]    6 200 0.5 0.5 0.1 0.1 1.0
 [7,]    7 200 0.3 0.3 0.1 0.1 0.0
 [8,]    8 200 0.7 0.7 0.1 0.1 0.0
 [9,]    9 200 0.3 0.5 0.1 0.1 0.0
[10,]   10 200 0.5 0.3 0.1 0.1 0.0
[11,]   11 200 0.4 0.6 0.1 0.2 1.0
[12,]   12 200 0.6 0.4 0.1 0.2 1.0
> 
> 
> 
> cleanEx(); ..nameEx <- "luo.ld.power"
> 
> ### * luo.ld.power
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: luo.ld.power
> ### Title: Classical deterministic power calculation for association
> ###   studies to detect linkage disequilibrium
> ### Aliases: luo.ld.power
> ### Keywords: design
> 
> ### ** Examples
> 
> data(luo.ld.populations)
> options(digits=3)
> powers <- numeric(nrow(luo.ld.populations))
> for(ii in 1:nrow(luo.ld.populations)){
+    cat("ii=",ii,"\n")
+    powers[ii] <- luo.ld.power(n=luo.ld.populations[ii,"n"],
+                               p=luo.ld.populations[ii,"p"],
+                               q=luo.ld.populations[ii,"q"],
+                               D=luo.ld.populations[ii,"D"],
+                               h2=luo.ld.populations[ii,"h2"],
+                               phi=luo.ld.populations[ii,"phi"],
+                               Vp=100,
+                               alpha=0.05)
+ }
ii= 1 
ems.beta = 178 , ems.omega = 98.4 , power = 0.181 
ii= 2 
ems.beta = 258 , ems.omega = 98.4 , power = 0.337 
ii= 3 
ems.beta = 730 , ems.omega = 93.6 , power = 0.915 
ii= 4 
ems.beta = 415 , ems.omega = 96.8 , power = 0.616 
ii= 5 
ems.beta = 243 , ems.omega = 98.5 , power = 0.308 
ii= 6 
ems.beta = 213 , ems.omega = 98.8 , power = 0.251 
ii= 7 
ems.beta = 325 , ems.omega = 97.7 , power = 0.464 
ii= 8 
ems.beta = 325 , ems.omega = 97.7 , power = 0.464 
ii= 9 
ems.beta = 287 , ems.omega = 98.1 , power = 0.393 
ii= 10 
ems.beta = 288 , ems.omega = 98.1 , power = 0.395 
ii= 11 
ems.beta = 318 , ems.omega = 97.8 , power = 0.452 
ii= 12 
ems.beta = 371 , ems.omega = 97.2 , power = 0.545 
> cbind(luo.ld.populations,power=powers)
      pop.   n   p   q   D  h2 phi power
 [1,]    1 100 0.5 0.5 0.1 0.1 0.0 0.181
 [2,]    2 200 0.5 0.5 0.1 0.1 0.0 0.337
 [3,]    3 200 0.5 0.5 0.2 0.1 0.0 0.915
 [4,]    4 200 0.5 0.5 0.1 0.2 0.0 0.616
 [5,]    5 200 0.5 0.5 0.1 0.1 0.5 0.308
 [6,]    6 200 0.5 0.5 0.1 0.1 1.0 0.251
 [7,]    7 200 0.3 0.3 0.1 0.1 0.0 0.464
 [8,]    8 200 0.7 0.7 0.1 0.1 0.0 0.464
 [9,]    9 200 0.3 0.5 0.1 0.1 0.0 0.393
[10,]   10 200 0.5 0.3 0.1 0.1 0.0 0.395
[11,]   11 200 0.4 0.6 0.1 0.2 1.0 0.452
[12,]   12 200 0.6 0.4 0.1 0.2 1.0 0.545
> 
> 
> 
> cleanEx(); ..nameEx <- "oneway.bf.alpha"
> 
> ### * oneway.bf.alpha
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: oneway.bf.alpha
> ### Title: Correspondence between significance levels and Bayes factors for
> ###   effects of marker genotype classes.
> ### Aliases: oneway.bf.alpha oneway.bf.alpha.required
> ### Keywords: htest models
> 
> ### ** Examples
> 
> # calculations for Table 4 in the manuscript
> data(luo.ld.populations)
> Bs <- numeric(nrow(luo.ld.populations))
> n.Bf20s <- numeric(nrow(luo.ld.populations))
> ns <- c(seq(200,400,by=25),,450,seq(500,4000,by=100))
> powers <- numeric(length(ns))
> alphas <- numeric(length(ns))
> P.Bf20s <- numeric(length(ns))
> for(ii in 1:nrow(luo.ld.populations)){
+    cat("ii=",ii,"\n")
+    powers[ii] <- luo.ld.power(n=luo.ld.populations[ii,"n"],
+                               p=luo.ld.populations[ii,"p"],
+                               q=luo.ld.populations[ii,"q"],
+                               D=luo.ld.populations[ii,"D"],
+                               h2=luo.ld.populations[ii,"h2"],
+                               phi=luo.ld.populations[ii,"phi"],
+                               Vp=100,
+                               alpha=0.05)
+ 
+    p1 <- luo.ld.populations[ii,"p"]
+    Bs[ii] <- oneway.bf.alpha(n=luo.ld.populations[ii,"n"],
+                  group.sizes=c(p1^2,2*p1*(1-p1),(1-p1)^2)*
+                                  luo.ld.populations[ii,"n"])
+    for(jj in seq(along=ns)){
+      alphas[jj] <- oneway.bf.alpha.required(ns[jj],
+                        group.sizes=c(p1^2,2*p1*(1-p1),(1-p1)^2)*ns[jj],Bf=20)
+      P.Bf20s[jj] <- luo.ld.power(n=ns[jj],
+                                  p=luo.ld.populations[ii,"p"],
+                                  luo.ld.populations[ii,"q"],
+                                  D=luo.ld.populations[ii,"D"],
+                                  h2=luo.ld.populations[ii,"h2"],
+                                  phi=luo.ld.populations[ii,"phi"],
+                                  Vp=100,
+                                  alpha=alphas[jj],
+                                  print.it=FALSE)
+    }
+    n.Bf20s[ii] <- approx(P.Bf20s,ns,xout=0.9)$y
+    cat("n =",n.Bf20s[ii],"\n")
+ }
ii= 1 
ems.beta = 178 , ems.omega = 98.4 , power = 0.181 
n = 1837 
ii= 2 
ems.beta = 258 , ems.omega = 98.4 , power = 0.337 
n = 1837 
ii= 3 
ems.beta = 730 , ems.omega = 93.6 , power = 0.915 
Warning in approx(P.Bf20s, ns, xout = 0.9) : 
	 collapsing to unique 'x' values
n = 381 
ii= 4 
ems.beta = 415 , ems.omega = 96.8 , power = 0.616 
n = 850 
ii= 5 
ems.beta = 243 , ems.omega = 98.5 , power = 0.308 
n = 2047 
ii= 6 
ems.beta = 213 , ems.omega = 98.8 , power = 0.251 
n = 2640 
ii= 7 
ems.beta = 325 , ems.omega = 97.7 , power = 0.464 
n = 1211 
ii= 8 
ems.beta = 325 , ems.omega = 97.7 , power = 0.464 
n = 1211 
ii= 9 
ems.beta = 287 , ems.omega = 98.1 , power = 0.393 
n = 1476 
ii= 10 
ems.beta = 288 , ems.omega = 98.1 , power = 0.395 
n = 1513 
ii= 11 
ems.beta = 318 , ems.omega = 97.8 , power = 0.452 
n = 1259 
ii= 12 
ems.beta = 371 , ems.omega = 97.2 , power = 0.545 
n = 995 
> cbind(luo.ld.populations,powers,n.Bf20s)
Warning in cbind(luo.ld.populations, powers, n.Bf20s) : 
	 number of rows of result
	is not a multiple of vector length (arg 2)
      pop.   n   p   q   D  h2 phi powers n.Bf20s
 [1,]    1 100 0.5 0.5 0.1 0.1 0.0  0.181    1837
 [2,]    2 200 0.5 0.5 0.1 0.1 0.0  0.337    1837
 [3,]    3 200 0.5 0.5 0.2 0.1 0.0  0.915     381
 [4,]    4 200 0.5 0.5 0.1 0.2 0.0  0.616     850
 [5,]    5 200 0.5 0.5 0.1 0.1 0.5  0.308    2047
 [6,]    6 200 0.5 0.5 0.1 0.1 1.0  0.251    2640
 [7,]    7 200 0.3 0.3 0.1 0.1 0.0  0.464    1211
 [8,]    8 200 0.7 0.7 0.1 0.1 0.0  0.464    1211
 [9,]    9 200 0.3 0.5 0.1 0.1 0.0  0.393    1476
[10,]   10 200 0.5 0.3 0.1 0.1 0.0  0.395    1513
[11,]   11 200 0.4 0.6 0.1 0.2 1.0  0.452    1259
[12,]   12 200 0.6 0.4 0.1 0.2 1.0  0.545     995
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  29.0 3.53 34.7 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')