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("hapassoc-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('hapassoc')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "PreHap"
> 
> ### * PreHap
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pre.hapassoc
> ### Title: Pre-process the data before fitting it with hapassoc
> ### Aliases: pre.hapassoc
> ### Keywords: methods
> 
> ### ** Examples
> 
> #First example data set has single-locus genotypes in "allelic format"
> data(hypoDat)
> example.pre.hapassoc<-pre.hapassoc(hypoDat, numSNPs=3)
> 
> # To get the initial haplotype frequencies:
> example.pre.hapassoc$initFreq
      h000       h001       h010       h011       h100       h101       h110 
0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844 0.01081260 
      h111 
0.02148268 
> #      h000       h001       h010       h011       h100       h101       h110 
> #0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844 0.01081260 
> #      h111 
> #0.02148268 
> # The '001' haplotype is estimated to be the most frequent
> 
> example.pre.hapassoc$pooledHaplos
[1] "h101" "h110" "h111"
> # "h101" "h110" "h111"
> # These haplotypes are to be pooled in the design matrix for the risk model
> 
> names(example.pre.hapassoc$haploDM)
[1] "h000"   "h001"   "h010"   "h011"   "h100"   "pooled"
> # "h000"   "h001"   "h010"   "h011"   "h100"   "pooled"
> 
> ####
> #Second example data set has single-locus genotypes in "genotypic format"
> data(hypoDatGeno)
> example2.pre.hapassoc<-pre.hapassoc(hypoDatGeno, numSNPs=3, allelic=FALSE)
> 
> # To get the initial haplotype frequencies:
> example2.pre.hapassoc$initFreq
      hAAA       hAAC       hACA       hACC       hCAA       hCAC       hCCA 
0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844 0.01081260 
      hCCC 
0.02148268 
> #      hAAA       hAAC       hACA       hACC       hCAA       hCAC
> #0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844
> #      hCCA           hCCC 
> #0.01081260     0.02148268 
> # The 'hAAC' haplotype is estimated to be the most frequent
> 
> example2.pre.hapassoc$pooledHaplos
[1] "hCAC" "hCCA" "hCCC"
> #  "hCAC" "hCCA" "hCCC"
> # These haplotypes are to be pooled in the design matrix for the risk model
> 
> names(example2.pre.hapassoc$haploDM)
[1] "hAAA"   "hAAC"   "hACA"   "hACC"   "hCAA"   "pooled"
> #  "hAAA"   "hAAC"   "hACA"   "hACC"   "hCAA"   "pooled"
> 
> 
> 
> cleanEx(); ..nameEx <- "hapassoc"
> 
> ### * hapassoc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hapassoc
> ### Title: EM algorithm to fit maximum likelihood estimates of trait
> ###   associations with SNP haplotypes
> ### Aliases: hapassoc
> ### Keywords: methods
> 
> ### ** Examples
> 
> data(hypoDat)
> example.pre.hapassoc<-pre.hapassoc(hypoDat, 3)
> 
> example.pre.hapassoc$initFreq # look at initial haplotype frequencies
      h000       h001       h010       h011       h100       h101       h110 
0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844 0.01081260 
      h111 
0.02148268 
> #      h000       h001       h010       h011       h100       h101       h110 
> #0.25179111 0.26050418 0.23606001 0.09164470 0.10133627 0.02636844 0.01081260 
> #      h111 
> #0.02148268 
> 
> names(example.pre.hapassoc$haploDM)
[1] "h000"   "h001"   "h010"   "h011"   "h100"   "pooled"
> # "h000"   "h001"   "h010"   "h011"   "h100"   "pooled"
> 
> # Columns of the matrix haploDM score the number of copies of each haplotype 
> # for each pseudo-individual.
> 
> # Logistic regression for a multiplicative odds model having as the baseline 
> # group homozygotes '001/001' for the most common haplotype
> 
> example.regr <- hapassoc(affected ~ attr + h000+ h010 + h011 + h100 + pooled,
+                   example.pre.hapassoc, family=binomial())
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
> 
> # Logistic regression with separate effects for 000 homozygotes, 001 homozygotes 
> # and 000/001 heterozygotes
> 
> example2.regr <- hapassoc(affected ~ attr + I(h000==2) + I(h001==2) +
+                    I(h000==1 & h001==1), example.pre.hapassoc, family=binomial())
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "summaryHap"
> 
> ### * summaryHap
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.hapassoc
> ### Title: Summarize results of the hapassoc function
> ### Aliases: summary.hapassoc
> ### Keywords: methods
> 
> ### ** Examples
> 
> data(hypoDat)
> example.pre.hapassoc<-pre.hapassoc(hypoDat, 3)
> example.regr <- hapassoc(affected ~ attr + h000+ h010 + h011 + h100 + pooled,
+                      example.pre.hapassoc, family=binomial())
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Warning in eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
> 
> # Summarize the results:
> summary.hapassoc(example.regr) # or just summary(example.regr)
$coefficients
               Estimate Std. Error      zscore   Pr(>|z|)
(Intercept) -1.24304712  0.7826121 -1.58833114 0.11221148
attr         0.74117082  0.2919206  2.53894658 0.01111868
h000         1.15178187  0.5942128  1.93833229 0.05258270
h010        -0.59528765  0.6570407 -0.90601329 0.36492882
h011        -0.03426242  0.9155250 -0.03742379 0.97014710
h100        -0.85628090  1.0199383 -0.83954187 0.40116530
pooled       0.38745841  0.8775181  0.44153893 0.65882289

$frequencies
         Estimate Std. Error
f.h000 0.26732097 0.03930403
f.h001 0.25180046 0.03865720
f.h010 0.21981267 0.03877367
f.h011 0.10106590 0.02949351
f.h100 0.09501225 0.02369630
f.h101 0.02586632 0.01412337
f.h110 0.01785410 0.01384235
f.h111 0.02126733 0.01247509

$dispersion
[1] 1

> 
> # Results:
> #$coefficients
> #               Estimate Std. Error      zscore   Pr(>|z|)
> #(Intercept) -1.24114270  0.7820977 -1.58694079 0.11252606
> #attr         0.74036920  0.2918205  2.53707057 0.01117844
> #h000         1.14968352  0.5942542  1.93466627 0.05303126
> #h010        -0.59318434  0.6569672 -0.90291311 0.36657201
> #h011        -0.03615243  0.9161959 -0.03945928 0.96852422
> #h100        -0.85329292  1.0203105 -0.83630709 0.40298217
> #pooled       0.38516864  0.8784283  0.43847478 0.66104215
> #
> #$frequencies
> #         Estimate Std. Error
> #f.h000 0.26716394 0.03933158
> #f.h001 0.25191674 0.03866739
> #f.h010 0.21997138 0.03881578
> #f.h011 0.10094795 0.02949617
> #f.h100 0.09507014 0.02371878
> #f.h101 0.02584918 0.01411881
> #f.h110 0.01779455 0.01386080
> #f.h111 0.02128613 0.01247265
> #
> #$dispersion
> #[1] 1
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  4.57 0.08 5.5 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')