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