R version 4.4.1 (2024-06-14) -- "Race for Your Life" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu 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. Natural language support but running in an English locale 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 an HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored] > > ## package: haplo.stats > ## test script: haplo.glm > ## created: 11/23/2011 > > ## settings > verbose=TRUE > require(haplo.stats) Loading required package: haplo.stats Loading required package: arsenal > Sys.setlocale("LC_COLLATE", "C") [1] "C" > Sys.getlocale('LC_COLLATE') [1] "C" > > > if(verbose) cat("setting up data...\n") setting up data... > > > # prepare the hla dataset, > # runs a lot longer, and MS alleles don't all start w/ 1, 2... > label <-c("DQB","DRB","B") > > data(hla.demo) > > y <- hla.demo$resp > y.bin <- 1*(hla.demo$resp.cat=="low") > geno <- as.matrix(hla.demo[,c(17,18,21:24)]) > geno <- setupGeno(geno, miss.val=c(0,NA)) > > # geno now has an attribute 'unique.alleles' which must be passed to > # haplo.glm as allele.lev=attributes(geno)$unique.alleles, see below > > hla.data <- data.frame(geno=geno, age=hla.demo$age, male=hla.demo$male, + y=y, y.bin=y.bin) > > seed <- c(17, 53, 1, 40, 37, 0, 62, 56, 5, 52, 12, 1) > > > if(verbose) cat("fit a binary trait\n") fit a binary trait > set.seed(seed) > fit.hla.bin <- haplo.glm(y.bin ~ male + geno, family = binomial, + na.action="na.geno.keep", data=hla.data, locus.label=label, + control = haplo.glm.control(haplo.min.count=8)) > > y.bin <- 1*(hla.demo$resp.cat=="low") > y.bin[2] <- NA > geno.hla <- as.matrix(hla.demo[,c(17,18,21:24)]) > geno.hla[2,5] <- 2 > geno.hla[3,] <- rep(NA, 6) > geno.hla <- setupGeno(geno.hla, miss.val=c(0,NA)) > > my.hla <- data.frame(geno.hla=geno.hla, age=hla.demo$age, male=hla.demo$male, + y=y, y.bin=y.bin) > > if(verbose) cat(" hla binary trait with subject that are removed\n") hla binary trait with subject that are removed > > set.seed(seed) > fit.hla.miss <- haplo.glm(y.bin ~ male + geno.hla, family = binomial, + na.action="na.geno.keep", + data=my.hla, locus.label=label, + control = haplo.glm.control(haplo.min.count=8)) > > > if(verbose) cat(" gaussian with covariates, additive\n") gaussian with covariates, additive > > set.seed(seed) > fit.hla.gaus.gender <- haplo.glm(y ~ male + geno, family = gaussian, + na.action="na.geno.keep", + data=hla.data, locus.label=label, + control = haplo.glm.control(haplo.min.count=5)) > > > if(verbose) cat("SNAP data with resp and resp with added variance\n") SNAP data with resp and resp with added variance > snapDF <- read.table("snapData.csv",header=TRUE, sep=",", stringsAsFactors=FALSE) > > geno.rec <- setupGeno(snapDF[,-c(1:9)]) > snap.data <- data.frame(resp=hla.demo$resp, respvar=hla.demo$resp*100, geno=geno.rec) > > set.seed(seed) > fit.resp.hla <- haplo.glm(resp~geno, trait.type="gaussian",data=snap.data) > > set.seed(seed) > fit.respvar.hla <- haplo.glm(respvar~geno, trait.type="gaussian",data=snap.data) > > > cat("summary function\n") summary function > > print(summary(fit.hla.bin),digits=3) Call: haplo.glm(formula = y.bin ~ male + geno, family = binomial, data = hla.data, na.action = "na.geno.keep", locus.label = label, control = haplo.glm.control(haplo.min.count = 8)) Coefficients: coef se t.stat pval (Intercept) 1.546 0.655 2.361 0.02 male -0.480 0.331 -1.452 0.15 geno.17 -0.723 0.801 -0.902 0.37 geno.34 0.364 0.680 0.536 0.59 geno.77 -0.988 0.733 -1.349 0.18 geno.78 -1.409 0.854 -1.650 0.10 geno.100 -2.591 1.128 -2.297 0.02 geno.138 -2.716 0.852 -3.186 0.00 geno.rare -1.261 0.354 -3.565 0.00 (Dispersion parameter for binomial family taken to be 1) Null deviance: 263.50 on 219 degrees of freedom Residual deviance: 233.46 on 211 degrees of freedom AIC: 251.1 Number of Fisher Scoring iterations: 61 Haplotypes: DQB DRB B hap.freq geno.17 21 7 44 0.0230 geno.34 31 4 44 0.0284 geno.77 32 4 60 0.0306 geno.78 32 4 62 0.0235 geno.100 51 1 35 0.0298 geno.138 62 2 7 0.0518 geno.rare * * * 0.7088 haplo.base 21 3 8 0.1041 > > > cat("fitted values for hlabin, hla-gaussian, hla-gaussian-hi-variance\n") fitted values for hlabin, hla-gaussian, hla-gaussian-hi-variance > > print(fitted(fit.hla.bin)[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 0.1673 0.1460 0.1890 0.3513 0.2736 0.5849 0.0808 0.1890 0.1383 0.0581 0.7019 12 13 14 15 16 17 18 19 20 0.2736 0.2736 0.0719 0.2736 0.0581 0.2736 0.0808 0.0516 0.1890 > > print(fitted(fit.resp.hla)[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 2.10 2.10 1.58 1.84 1.72 1.85 1.83 1.83 2.09 2.09 1.83 2.09 1.83 1.83 1.84 2.09 17 18 19 20 1.87 1.83 1.88 1.85 > print(fitted(fit.respvar.hla)[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 210 210 158 184 172 185 183 183 209 209 183 209 183 183 184 209 187 183 188 185 > > cat("vcov for hlabin, hla-gaussian, hla-gaussian-hi-variance\n") vcov for hlabin, hla-gaussian, hla-gaussian-hi-variance > > print(vcov(fit.hla.bin)[1:20,1:20],digits=3) (Intercept) male geno.17 geno.34 geno.77 geno.78 (Intercept) 4.29e-01 -7.29e-02 -2.37e-01 -1.80e-01 -1.83e-01 -1.83e-01 male -7.29e-02 1.09e-01 2.76e-02 -1.13e-02 -1.74e-02 1.52e-02 geno.17 -2.37e-01 2.76e-02 6.42e-01 1.08e-01 6.36e-02 1.12e-01 geno.34 -1.80e-01 -1.13e-02 1.08e-01 4.62e-01 7.42e-02 9.14e-02 geno.77 -1.83e-01 -1.74e-02 6.36e-02 7.42e-02 5.37e-01 8.95e-02 geno.78 -1.83e-01 1.52e-02 1.12e-01 9.14e-02 8.95e-02 7.30e-01 geno.100 -2.48e-01 4.86e-02 1.48e-01 1.11e-01 5.85e-02 1.14e-01 geno.138 -2.35e-01 3.14e-02 1.06e-01 6.04e-02 9.92e-02 8.99e-02 geno.rare -2.15e-01 1.20e-02 1.14e-01 9.40e-02 9.75e-02 8.48e-02 hap.17 6.51e-06 9.19e-07 -3.91e-06 -3.65e-06 -3.83e-06 -3.19e-06 hap.34 -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.77 -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.78 -6.82e-07 -3.28e-08 3.01e-07 2.95e-07 3.15e-07 2.24e-07 hap.100 5.14e-06 1.41e-06 6.12e-06 -3.08e-06 -4.34e-06 -2.13e-06 hap.138 1.40e-05 2.52e-06 -1.50e-05 -1.09e-05 -1.25e-05 -8.45e-06 hap.rare -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.rare -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.rare -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.rare -4.06e-08 3.31e-08 1.09e-08 -1.24e-09 -2.57e-09 8.52e-10 hap.rare -8.12e-08 6.62e-08 2.19e-08 -2.48e-09 -5.14e-09 1.70e-09 geno.100 geno.138 geno.rare hap.17 hap.34 hap.77 (Intercept) -2.48e-01 -2.35e-01 -2.15e-01 6.51e-06 -4.06e-08 -4.06e-08 male 4.86e-02 3.14e-02 1.20e-02 9.19e-07 3.31e-08 3.31e-08 geno.17 1.48e-01 1.06e-01 1.14e-01 -3.91e-06 1.09e-08 1.09e-08 geno.34 1.11e-01 6.04e-02 9.40e-02 -3.65e-06 -1.24e-09 -1.24e-09 geno.77 5.85e-02 9.92e-02 9.75e-02 -3.83e-06 -2.57e-09 -2.57e-09 geno.78 1.14e-01 8.99e-02 8.48e-02 -3.19e-06 8.52e-10 8.52e-10 geno.100 1.27e+00 1.16e-01 1.14e-01 -3.69e-06 1.37e-08 1.37e-08 geno.138 1.16e-01 7.27e-01 1.16e-01 -3.79e-06 8.44e-09 8.44e-09 geno.rare 1.14e-01 1.16e-01 1.25e-01 -3.69e-06 2.65e-08 2.65e-08 hap.17 -3.69e-06 -3.79e-06 -3.69e-06 5.52e-06 -1.21e-08 -1.21e-08 hap.34 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 5.15e-06 -1.17e-08 hap.77 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 -1.17e-08 5.15e-06 hap.78 3.15e-07 3.46e-07 4.74e-07 -1.22e-08 -1.17e-08 -1.17e-08 hap.100 8.13e-07 -3.34e-06 -3.58e-06 -3.04e-08 -2.96e-08 -2.96e-08 hap.138 -9.20e-06 -8.24e-06 -6.49e-06 -1.89e-08 -1.93e-08 -1.93e-08 hap.rare 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 -1.17e-08 -1.17e-08 hap.rare 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 -1.17e-08 -1.17e-08 hap.rare 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 -1.17e-08 -1.17e-08 hap.rare 1.37e-08 8.44e-09 2.65e-08 -1.21e-08 -1.17e-08 -1.17e-08 hap.rare 2.75e-08 1.69e-08 5.30e-08 -2.42e-08 -2.34e-08 -2.34e-08 hap.78 hap.100 hap.138 hap.rare hap.rare hap.rare (Intercept) -6.82e-07 5.14e-06 1.40e-05 -4.06e-08 -4.06e-08 -4.06e-08 male -3.28e-08 1.41e-06 2.52e-06 3.31e-08 3.31e-08 3.31e-08 geno.17 3.01e-07 6.12e-06 -1.50e-05 1.09e-08 1.09e-08 1.09e-08 geno.34 2.95e-07 -3.08e-06 -1.09e-05 -1.24e-09 -1.24e-09 -1.24e-09 geno.77 3.15e-07 -4.34e-06 -1.25e-05 -2.57e-09 -2.57e-09 -2.57e-09 geno.78 2.24e-07 -2.13e-06 -8.45e-06 8.52e-10 8.52e-10 8.52e-10 geno.100 3.15e-07 8.13e-07 -9.20e-06 1.37e-08 1.37e-08 1.37e-08 geno.138 3.46e-07 -3.34e-06 -8.24e-06 8.44e-09 8.44e-09 8.44e-09 geno.rare 4.74e-07 -3.58e-06 -6.49e-06 2.65e-08 2.65e-08 2.65e-08 hap.17 -1.22e-08 -3.04e-08 -1.89e-08 -1.21e-08 -1.21e-08 -1.21e-08 hap.34 -1.17e-08 -2.96e-08 -1.93e-08 -1.17e-08 -1.17e-08 -1.17e-08 hap.77 -1.17e-08 -2.96e-08 -1.93e-08 -1.17e-08 -1.17e-08 -1.17e-08 hap.78 5.19e-06 -2.97e-08 -1.93e-08 -1.17e-08 -1.17e-08 -1.17e-08 hap.100 -2.97e-08 1.53e-05 -2.29e-06 -2.96e-08 -2.96e-08 -2.96e-08 hap.138 -1.93e-08 -2.29e-06 1.17e-05 -1.93e-08 -1.93e-08 -1.93e-08 hap.rare -1.17e-08 -2.96e-08 -1.93e-08 5.15e-06 -1.17e-08 -1.17e-08 hap.rare -1.17e-08 -2.96e-08 -1.93e-08 -1.17e-08 5.15e-06 -1.17e-08 hap.rare -1.17e-08 -2.96e-08 -1.93e-08 -1.17e-08 -1.17e-08 5.15e-06 hap.rare -1.17e-08 -2.96e-08 -1.93e-08 -1.17e-08 -1.17e-08 -1.17e-08 hap.rare -2.35e-08 -5.92e-08 -3.87e-08 -2.34e-08 -2.34e-08 -2.34e-08 hap.rare hap.rare (Intercept) -4.06e-08 -8.12e-08 male 3.31e-08 6.62e-08 geno.17 1.09e-08 2.19e-08 geno.34 -1.24e-09 -2.48e-09 geno.77 -2.57e-09 -5.14e-09 geno.78 8.52e-10 1.70e-09 geno.100 1.37e-08 2.75e-08 geno.138 8.44e-09 1.69e-08 geno.rare 2.65e-08 5.30e-08 hap.17 -1.21e-08 -2.42e-08 hap.34 -1.17e-08 -2.34e-08 hap.77 -1.17e-08 -2.34e-08 hap.78 -1.17e-08 -2.35e-08 hap.100 -2.96e-08 -5.92e-08 hap.138 -1.93e-08 -3.87e-08 hap.rare -1.17e-08 -2.34e-08 hap.rare -1.17e-08 -2.34e-08 hap.rare -1.17e-08 -2.34e-08 hap.rare 5.15e-06 -2.34e-08 hap.rare -2.34e-08 1.03e-05 > > print(vcov(fit.resp.hla),digits=3) (Intercept) geno.1 geno.4 geno.6 geno.7 hap.1 (Intercept) 2.84e-02 -2.07e-02 -1.50e-02 -1.41e-02 -2.01e-02 1.05e-04 geno.1 -2.07e-02 4.04e-02 9.22e-03 5.81e-03 1.41e-02 -3.19e-05 geno.4 -1.50e-02 9.22e-03 2.34e-02 3.46e-03 8.10e-03 -9.00e-05 geno.6 -1.41e-02 5.81e-03 3.46e-03 4.65e-02 5.10e-03 3.17e-04 geno.7 -2.01e-02 1.41e-02 8.10e-03 5.10e-03 2.87e-02 -5.58e-05 hap.1 1.05e-04 -3.19e-05 -9.00e-05 3.17e-04 -5.58e-05 3.26e-04 hap.4 -2.45e-06 2.66e-05 -2.77e-05 1.35e-05 -3.90e-07 -5.37e-05 hap.6 -9.37e-05 -8.67e-06 1.17e-04 -3.83e-04 6.44e-05 -1.04e-04 hap.7 9.37e-05 8.67e-06 -1.17e-04 3.83e-04 -6.44e-05 1.29e-05 hap.4 hap.6 hap.7 (Intercept) -2.45e-06 -9.37e-05 9.37e-05 geno.1 2.66e-05 -8.67e-06 8.67e-06 geno.4 -2.77e-05 1.17e-04 -1.17e-04 geno.6 1.35e-05 -3.83e-04 3.83e-04 geno.7 -3.90e-07 6.44e-05 -6.44e-05 hap.1 -5.37e-05 -1.04e-04 1.29e-05 hap.4 3.36e-04 -4.21e-05 -8.82e-05 hap.6 -4.21e-05 2.90e-04 -1.30e-04 hap.7 -8.82e-05 -1.30e-04 4.68e-04 > print(vcov(fit.respvar.hla),digits=3) (Intercept) geno.1 geno.4 geno.6 geno.7 hap.1 (Intercept) 2.84e+02 -2.07e+02 -1.50e+02 -1.41e+02 -2.01e+02 1.05e-02 geno.1 -2.07e+02 4.04e+02 9.22e+01 5.81e+01 1.41e+02 -3.19e-03 geno.4 -1.50e+02 9.22e+01 2.34e+02 3.46e+01 8.10e+01 -9.00e-03 geno.6 -1.41e+02 5.81e+01 3.46e+01 4.65e+02 5.10e+01 3.17e-02 geno.7 -2.01e+02 1.41e+02 8.10e+01 5.10e+01 2.87e+02 -5.58e-03 hap.1 1.05e-02 -3.19e-03 -9.00e-03 3.17e-02 -5.58e-03 3.26e-04 hap.4 -2.45e-04 2.66e-03 -2.77e-03 1.35e-03 -3.90e-05 -5.37e-05 hap.6 -9.37e-03 -8.67e-04 1.17e-02 -3.83e-02 6.44e-03 -1.04e-04 hap.7 9.37e-03 8.67e-04 -1.17e-02 3.83e-02 -6.44e-03 1.29e-05 hap.4 hap.6 hap.7 (Intercept) -2.45e-04 -9.37e-03 9.37e-03 geno.1 2.66e-03 -8.67e-04 8.67e-04 geno.4 -2.77e-03 1.17e-02 -1.17e-02 geno.6 1.35e-03 -3.83e-02 3.83e-02 geno.7 -3.90e-05 6.44e-03 -6.44e-03 hap.1 -5.37e-05 -1.04e-04 1.29e-05 hap.4 3.36e-04 -4.21e-05 -8.82e-05 hap.6 -4.21e-05 2.90e-04 -1.30e-04 hap.7 -8.82e-05 -1.30e-04 4.68e-04 > > cat("residuals for hlabin, hla-gaussian, hla-gaussian-hi-variance\n") residuals for hlabin, hla-gaussian, hla-gaussian-hi-variance > > print(residuals(fit.hla.bin, type="deviance")[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 -0.605 -0.562 -0.647 -0.930 1.610 -1.326 -0.411 -0.647 -0.546 -0.346 -1.556 12 13 14 15 16 17 18 19 20 -0.800 -0.800 -0.386 -0.800 -0.346 -0.800 -0.411 -0.326 -0.647 > print(residuals(fit.hla.bin, type="pearson")[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 -0.448 -0.414 -0.483 -0.736 1.629 -1.187 -0.297 -0.483 -0.401 -0.248 -1.535 12 13 14 15 16 17 18 19 20 -0.614 -0.614 -0.278 -0.614 -0.248 -0.614 -0.297 -0.233 -0.483 > print(residuals(fit.hla.bin, type="response")[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 -0.1673 -0.1460 -0.1890 -0.3513 0.7264 -0.5849 -0.0808 -0.1890 -0.1383 -0.0581 11 12 13 14 15 16 17 18 19 20 -0.7019 -0.2736 -0.2736 -0.0719 -0.2736 -0.0581 -0.2736 -0.0808 -0.0516 -0.1890 > > print(residuals(fit.resp.hla)[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 1.5972 -0.9248 0.0778 0.1307 -0.9213 -0.6994 -0.5256 0.2524 -0.9300 -0.5340 11 12 13 14 15 16 17 18 19 20 -0.6046 0.2880 0.5964 0.3824 0.5287 0.8420 0.6858 -0.7366 0.7558 -0.1524 > print(residuals(fit.respvar.hla)[1:20],digits=3) 1 2 3 4 5 6 7 8 9 10 11 159.72 -92.48 7.78 13.07 -92.13 -69.94 -52.56 25.24 -93.00 -53.40 -60.46 12 13 14 15 16 17 18 19 20 28.80 59.64 38.24 52.87 84.20 68.58 -73.66 75.58 -15.24 > > > proc.time() user system elapsed 14.182 0.351 15.251