[Genabel-commits] r1105 - in branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests: . test_cox_ngpreds2
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 14 17:46:22 CET 2013
Author: lckarssen
Date: 2013-02-14 17:46:21 +0100 (Thu, 14 Feb 2013)
New Revision: 1105
Added:
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_pacox.R
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_palinear.R
Log:
Added and updated Yurii's R checks from the pacoxfix v0.1.3 branch
Added: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_pacox.R
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_pacox.R (rev 0)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_pacox.R 2013-02-14 16:46:21 UTC (rev 1105)
@@ -0,0 +1,111 @@
+library(survival)
+
+###
+# load the data
+###
+
+# load phenotypic data
+pheno <- read.table("../../examples/coxph_data.txt", head=TRUE,
+ string=FALSE)
+
+# load genetic DOSE data
+dose <- read.table("../../examples/test.mldose", head=FALSE,
+ string=FALSE)
+# remove "1->" from the names of dose-IDs
+idNames <- dose[, 1]
+idNames <- sub("[0-9]+->", "", idNames)
+dose[,1] <- idNames
+# check consistency of names
+table(dose[, 1] == pheno[, 1])
+
+# load genetic PROB data
+prob <- read.table("../../examples/test.mlprob", head=FALSE, string=FALSE)
+# remove "1->" from the names of prob-IDs
+idNames <- prob[, 1]
+idNames <- sub("[0-9]+->", "", idNames)
+prob[,1] <- idNames
+# check consistency of names
+table(prob[,1] == pheno[,1])
+
+# check consistency DOSE <-> PROB
+doseFromProb <- matrix(NA, ncol=dim(dose)[2], nrow=dim(dose)[1])
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ doseFromProb[, i] <- prob[, indexHom] * 2 + prob[, indexHet]
+ print( table( ( dose[,i] - doseFromProb[,i] ) > 1e-8 ) )
+}
+
+###
+# run analysis
+###
+
+# run ProbABEL
+system("cd ../../examples/; sh example_cox.sh; cd -")
+resPaAddDose <-
+ read.table("../../examples/coxph_dose_add.out.txt",
+ head=TRUE)[, c("beta_SNP_add","sebeta_SNP_add")]
+resPaAddProb <-
+ read.table("../../examples/coxph_prob_add.out.txt",
+ head=TRUE)[, c("beta_SNP_addA1","sebeta_SNP_addA1")]
+resPaDom <-
+ read.table("../../examples/coxph_prob_domin.out.txt",
+ head=TRUE)[, c("beta_SNP_domA1","sebeta_SNP_domA1")]
+resPaRec <-
+ read.table("../../examples/coxph_prob_recess.out.txt",
+ head=TRUE)[, c("beta_SNP_recA1","sebeta_SNP_recA1")]
+resPaOdom <-
+ read.table("../../examples/coxph_prob_over_domin.out.txt",
+ head=TRUE)[, c("beta_SNP_odom","sebeta_SNP_odom")]
+
+attach(pheno)
+
+# run analysis on dose
+emptyRes <- resPaAddDose
+emptyRes[] <- NA
+resRAddDose <- resRAddProb <- resRDom <- resRRec <- resROdom <- emptyRes
+for (i in 3:dim(dose)[2]) {
+ smr <- summary( coxph(Surv(fupt_chd, chd) ~ sex + age +
+ othercov + dose[,i] ) )
+ resRAddDose[i-2,] <- smr$coef[4, c(1,3)]
+ smr <- summary( coxph(Surv(fupt_chd, chd) ~ sex + age +
+ othercov + doseFromProb[,i] ) )
+ resRAddProb[i-2,] <- smr$coef[4, c(1,3)]
+}
+cat("additive model (dose):\n")
+resPaAddDose/resRAddDose
+cat("additive model (prob):\n")
+resPaAddProb/resRAddProb
+
+cat("dominant model (prob):\n")
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHom] + prob[,indexHet]
+ smr <- summary( coxph(Surv(fupt_chd, chd) ~ sex + age +
+ othercov + regProb ) )
+ resRDom[i-2,] <- smr$coef[4, c(1,3)]
+}
+resPaDom/resRDom
+
+cat("recessive model (prob):\n")
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHom]
+ smr <- summary( coxph(Surv(fupt_chd, chd) ~ sex + age +
+ othercov + regProb ) )
+ resRRec[i-2,] <- smr$coef[4, c(1,3)]
+}
+resPaRec/resRRec
+
+cat("over-dominant model (prob):\n")
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[,indexHet]
+ smr <- summary( coxph(Surv(fupt_chd, chd) ~ sex + age +
+ othercov + regProb ) )
+ resROdom[i-2,] <- smr$coef[4, c(1,3)]
+}
+resPaOdom/resROdom
Added: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_palinear.R
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_palinear.R (rev 0)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/tests/test_cox_ngpreds2/run_models_in_R_palinear.R 2013-02-14 16:46:21 UTC (rev 1105)
@@ -0,0 +1,74 @@
+###
+# load the data
+###
+
+# load phenotypic data
+pheno <- read.table("../../examples/height.txt", head=TRUE, string=FALSE)
+
+# load genetic DOSE data
+dose <- read.table("../../examples/test.mldose.2", head=FALSE, string=FALSE)
+# remove "1->" from the names of dose-IDs
+idNames <- dose[, 1]
+idNames <- sub("[0-9]+->", "", idNames)
+dose[, 1] <- idNames
+# check consistency of names
+table(dose[, 1] == pheno[, 1])
+
+# load genetic PROB data
+prob <- read.table("../../examples/test.mlprob", head=FALSE, string=FALSE)
+# remove "1->" from the names of prob-IDs
+idNames <- prob[, 1]
+idNames <- sub("[0-9]+->", "", idNames)
+prob[, 1] <- idNames
+# check consistency of names
+table(prob[, 1] == pheno[, 1])
+
+# check consistency DOSE <-> PROB
+doseFromProb <- matrix(NA, ncol=dim(dose)[2], nrow=dim(dose)[1])
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ doseFromProb[, i] <- prob[, indexHom] * 2 + prob[, indexHet]
+ print( table( ( dose[, i] - doseFromProb[, i] ) > 1e-8 ) )
+}
+
+###
+# run analysis
+###
+
+attach(pheno)
+
+# run analysis on dose
+for (i in 3:dim(dose)[2]) {
+ smr <- summary( lm( height ~ sex + age + dose[, i] ) )
+ print( smr$coef[4, 1:2] )
+ smr <- summary( lm( height ~ sex + age + doseFromProb[, i] ) )
+ print( smr$coef[4, 1:2] )
+}
+
+# dominant model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[, indexHom] + prob[, indexHet]
+ smr <- summary( lm( height ~ sex + age + regProb ) )
+ print( smr$coef[4, 1:2] )
+}
+
+# recessive model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[, indexHom]
+ smr <- summary( lm( height ~ sex + age + regProb ) )
+ print( smr$coef[4, 1:2] )
+}
+
+# over-dominant model
+for (i in 3:dim(dose)[2]) {
+ indexHom <- 3 + ( i - 3 ) * 2
+ indexHet <- indexHom + 1
+ regProb <- prob[, indexHet]
+ smr <- summary( lm( height ~ sex + age + regProb ) )
+ print( smr$coef[4, 1:2] )
+}
More information about the Genabel-commits
mailing list