[Genabel-commits] r1931 - in branches/ProbABEL-v0.4.4-hotfix/ProbABEL: checks checks/R-tests checks/inputfiles src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 17 14:54:45 CET 2015
Author: lckarssen
Date: 2015-03-17 14:54:45 +0100 (Tue, 17 Mar 2015)
New Revision: 1931
Added:
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh
Modified:
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/Makefile.am
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/Makefile.am
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_model_coxph.R
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/phedata.cpp
Log:
this commit fixes several bugs in ProbABEL:
- #1266: "pacoxph with no covariates"; I added a \mu of 1s to the X matrix. This doesn't affect the estimates for the SNPs and has the added benefit that the function coxph_data::coxph_data(phedata &phed, gendata &gend, const int snpnum) looks much more like it's cousin regdata::regdata(phedata &phed, gendata &gend, const int snpnum, const bool ext_is_interaction_excluded)).
- #6041: "ProbABEL's Cox module reports too many errors (beta may be infinite, setting beta and se to 'NaN') "; I added some more checks following the R Survial example code. I think the checks are now almost identical.
- This also contains a small "fix" in phedata.cpp in the way the model is printed on the screen. By moving the + to a different place we don't end up with things like height ~ mu + + snp when no covariates are present.
Implementation not related to the two bug fixes:
Adding a constant mean to the Cox PH regression I lead to "X matrix is singular" messages in ProbABEL and the R checks (luckily in both ;-)), Therefore, I had to adapt the R tests to disregard this message when related to the first variable. Similarly in the pacoxph code (in the estimate function).
Moreover, for the dominant model an error occurred in the R checks:
<simpleError in fitter(X, Y, strats, offset, init, control, weights = weights,
method = method, row.names(mf)): NA/NaN/Inf in foreign function call (arg 6)>
So I now catch errors in run_model_coxph.R as well and filter out the error I got.
I also added a copy of the text_cox.sh script, this time for the case without covariates.
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/Makefile.am
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/Makefile.am 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/Makefile.am 2015-03-17 13:54:45 UTC (rev 1931)
@@ -17,6 +17,7 @@
phenofiles = $(input_files_dir)/height.txt \
$(input_files_dir)/logist_data.txt \
$(input_files_dir)/coxph_data.txt \
+ $(input_files_dir)/coxph_data_nocovar.txt \
$(input_files_dir)/mmscore_pheno.PHE
genofiles = $(input_files_dir)/test.mldose \
@@ -49,7 +50,7 @@
check_SCRIPTS += test_bt.sh
endif
if BUILD_pacoxph
-check_SCRIPTS += test_cox.sh
+check_SCRIPTS += test_cox.sh test_cox_nocovar.sh
endif
AM_TESTS_ENVIRONMENT= \
@@ -125,6 +126,18 @@
coxph_prob_recess.out.txt coxph_prob_add.out.txt \
coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt \
+coxph_dose_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_over_domin.out.txt \
+coxph_prob_nocovar_fv_domin.out.txt \
+coxph_prob_nocovar_over_domin.out.txt \
+coxph_prob_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_recess.out.txt \
+coxph_prob_nocovar_domin.out.txt \
+coxph_prob_nocovar_recess.out.txt \
+coxph_prob_nocovar_add.out.txt \
+coxph_prob_nocovar_2df.out.txt \
+coxph_prob_nocovar_fv_2df.out.txt
dose_files = chr1.chunk1.dose chr1.chunk2.dose chr2.chunk1.dose \
chr2.chunk2.dose chr1.dose chr2.dose
@@ -147,4 +160,5 @@
cleanfiles_dose_check = test.mldose
CLEANFILES = $(cleanfiles_probabel_check) $(cleanfiles_dose_check) \
- $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox)
+ $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox) \
+ $(cleanfiles_cox_nocovar)
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/Makefile.am
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/Makefile.am 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/Makefile.am 2015-03-17 13:54:45 UTC (rev 1931)
@@ -104,5 +104,18 @@
coxph_prob_recess.out.txt coxph_prob_add.out.txt \
coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt \
+coxph_dose_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_over_domin.out.txt \
+coxph_prob_nocovar_fv_domin.out.txt \
+coxph_prob_nocovar_over_domin.out.txt \
+coxph_prob_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_recess.out.txt \
+coxph_prob_nocovar_domin.out.txt \
+coxph_prob_nocovar_recess.out.txt \
+coxph_prob_nocovar_add.out.txt \
+coxph_prob_nocovar_2df.out.txt \
+coxph_prob_nocovar_fv_2df.out.txt
+
CLEANFILES = $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) \
- $(cleanfiles_cox)
+ $(cleanfiles_cox) $(cleanfiles_cox_nocovar)
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_model_coxph.R
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_model_coxph.R 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_model_coxph.R 2015-03-17 13:54:45 UTC (rev 1931)
@@ -15,8 +15,8 @@
##' output can be compared to the ProbABEL output.
##' @author L.C. Larsen
run.model <- function(model0.txt, model.txt,
- snpcomponent1, snpcomponent2="1") {
-
+ snpcomponent1, snpcomponent2="1",
+ verbose=FALSE) {
if (snpcomponent2 != "1") {
## SNP component 2 is not constant: assume we run the 2df
## model.
@@ -28,15 +28,31 @@
resultR <- data.frame()
for (i in 3:dim(dose)[2]) {
+ if (verbose) print(paste("------- new iteration: i =", i))
indexHom <- 3 + ( i - 3 ) * 2
indexHet <- indexHom + 1
snp1 <- eval(parse(text=snpcomponent1))
snp2 <- eval(parse(text=snpcomponent2))
snp <- snp1 + snp2
+ mu <- rep(1., length(snp)) # Add a constant mean
noNA <- which( !is.na(snp) )
model.0 <- eval(parse(text=model0.txt))
+ ## Check the imputation R^2, if below threshold ProbABEL will
+ ## set the coefficients to NaN and skip the rest of this
+ ## iteration of the for loop.
+ rsq <- Rsq[i-2]
+ if (rsq < rsq.thresh ) {
+ if (twoDF) {
+ row <- c(rsq, NaN, NaN, NaN, NaN, NaN)
+ } else {
+ row <- c(rsq, NaN, NaN, NaN)
+ }
+ resultR <- rbind(resultR, row)
+ next
+ }
+
## Evaluate the model. The whole tryCatch is needed to catch
## problems with non-converging regression.
model = tryCatch({
@@ -49,11 +65,25 @@
eval(parse(text=model.txt)),
war)
)
- })
+ },
+ error = function(err) {
+ return(list(
+ "Can't calculate model; some error occurred",
+ err)
+ )
+ })
- if ( grepl("infinite", model[[2]]$message) |
- grepl("singular", model[[2]]$message) ) {
- ## The model did not converge, fill the coefficients with
+ if (verbose) print(model)
+
+ if ( grepl("Inf", model[[2]]$message) |
+ grepl("infinite", model[[2]]$message) |
+ grepl("iterations", model[[2]]$message) |
+ ( grepl("singular", model[[2]]$message) &
+ (regexpr("variable 1$", model[[2]]$message)[[1]] != 33)
+ )
+ ) {
+ ## The model did not converge or some other
+ ## errors/warnings occurred, fill the coefficients with
## NaNs
if (twoDF) {
smA1A2 <- c(NaN, NaN)
@@ -66,30 +96,21 @@
## No convergence problems, we can trust the
## coefficients.
coeff <- summary(model[[1]])$coefficients
+ # The SNP coefficient(s) is/are in the last row(s) of coeff
+ nr <- nrow(coeff)
if (twoDF) {
- smA1A2 <- coeff[4, c("coef", "se(coef)")]
- smA1A1 <- coeff[5, c("coef", "se(coef)")]
+ smA1A2 <- coeff[nr - 1, c("coef", "se(coef)")]
+ smA1A1 <- coeff[nr, c("coef", "se(coef)")]
} else {
- sm <- coeff[4, c("coef", "se(coef)")]
+ sm <- coeff[nr, c("coef", "se(coef)")]
}
lrt <- 2 * ( model[[1]]$loglik[2] - model.0$loglik[2] )
}
- ## Check the imputation R^2, if below threshold ProbABEL will
- ## set the coefficients to NaN.
- rsq <- Rsq[i-2]
if (twoDF) {
- if( rsq < rsq.thresh ) {
- row <- c(rsq, NaN, NaN, NaN, NaN, NaN)
- } else {
- row <- c(rsq, smA1A2[1], smA1A2[2], smA1A1[1], smA1A1[2], lrt)
- }
+ row <- c(rsq, smA1A2[1], smA1A2[2], smA1A1[1], smA1A1[2], lrt)
} else {
- if( rsq < rsq.thresh ) {
- row <- c(rsq, NaN, NaN, NaN)
- } else {
- row <- c(rsq, sm[1], sm[2], lrt)
- }
+ row <- c(rsq, sm[1], sm[2], lrt)
}
resultR <- rbind(resultR, row)
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2015-03-17 13:54:45 UTC (rev 1931)
@@ -53,8 +53,8 @@
source(paste0(srcdir, "run_model_coxph.R"))
model.fn.0 <-
- "coxph( Surv(fupt_chd, chd)[noNA] ~ sex[noNA] + age[noNA] + othercov[noNA] )"
-model.fn <- "coxph( Surv(fupt_chd, chd) ~ sex + age + othercov + snp1 )"
+ "coxph( Surv(fupt_chd, chd)[noNA] ~ mu[noNA] + sex[noNA] + age[noNA] + othercov[noNA] )"
+model.fn <- "coxph( Surv(fupt_chd, chd) ~ mu + sex + age + othercov + snp1 )"
## Additive model, dosages
snpdose <- "dose[, i]"
@@ -109,4 +109,85 @@
stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
cat("2df\n")
+
+
+cat("\nRun the same checks, but without any covariates (see bug #1266)\n")
+cat("Running ProbABEL...\t\t\t\t")
+tmp <- system(paste0("bash ", tests.path, "test_cox_nocovar.sh 2> /dev/null"),
+ intern=TRUE)
+cat("OK\n")
+
+dose.add.PA <- read.table("coxph_dose_nocovar_add.out.txt",
+ head=TRUE)[, colsAddDose]
+prob.add.PA <- read.table("coxph_prob_nocovar_add.out.txt",
+ head=TRUE)[, colsAddProb]
+prob.dom.PA <- read.table("coxph_prob_nocovar_domin.out.txt",
+ head=TRUE)[, colsDom]
+prob.rec.PA <- read.table("coxph_prob_nocovar_recess.out.txt",
+ head=TRUE)[, colsRec]
+prob.odom.PA <- read.table("coxph_prob_nocovar_over_domin.out.txt",
+ head=TRUE)[, colsOdom]
+prob.2df.PA <- read.table("coxph_prob_nocovar_2df.out.txt",
+ head=TRUE)[, cols2df]
+
+model.fn.0 <-
+ "coxph( Surv(fupt_chd, chd)[noNA] ~ mu[noNA] )"
+model.fn <- "coxph( Surv(fupt_chd, chd) ~ mu + snp1 )"
+
+
+cat("Comparing R output with ProbABEL output\t\t")
+
+## Additive model, dosages
+snpdose <- "dose[, i]"
+dose.add.R <- run.model(model.fn.0, model.fn, snpdose)
+colnames(dose.add.R) <- colsAddDose
+rownames(dose.add.R) <- NULL
+stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
+cat("additive ")
+
+
+## Additive model, probabilities
+snpprob <- "doseFromProb[, i]"
+prob.add.R <- run.model(model.fn.0, model.fn, snpprob)
+colnames(prob.add.R) <- colsAddProb
+rownames(prob.add.R) <- NULL
+stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
+cat("additive ")
+
+## dominant model
+snpprob <- "prob[, indexHom] + prob[, indexHet]"
+prob.dom.R <- run.model(model.fn.0, model.fn, snpprob)
+colnames(prob.dom.R) <- colsDom
+rownames(prob.dom.R) <- NULL
+stopifnot( all.equal(prob.dom.PA, prob.dom.R, tol=tol) )
+cat("dominant ")
+
+## recessive model
+snpprob <- "prob[, indexHom]"
+prob.rec.R <- run.model(model.fn.0, model.fn, snpprob)
+colnames(prob.rec.R) <- colsRec
+rownames(prob.rec.R) <- NULL
+stopifnot( all.equal(prob.rec.PA, prob.rec.R, tol=tol) )
+cat("recessive ")
+
+## over-dominant model
+snpprob <- "prob[, indexHet]"
+prob.odom.R <- run.model(model.fn.0, model.fn, snpprob)
+colnames(prob.odom.R) <- colsOdom
+rownames(prob.odom.R) <- NULL
+stopifnot( all.equal(prob.odom.PA, prob.odom.R, tol=tol) )
+cat("overdominant ")
+
+
+## 2df model
+model.fn <-
+ "coxph( Surv(fupt_chd, chd) ~ mu+ snp1 + snp2 )"
+snpd1 <- "prob[, indexHet]"
+snpd2 <- "prob[, indexHom]"
+prob.2df.R <- run.model(model.fn.0, model.fn, snpd1, snpd2)
+colnames(prob.2df.R) <- cols2df
+rownames(prob.2df.R) <- NULL
+stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
+cat("2df\n")
+
cat("\t\t\t\t\t\tOK\n")
Added: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt (rev 0)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt 2015-03-17 13:54:45 UTC (rev 1931)
@@ -0,0 +1,201 @@
+id fupt_chd chd
+id636728 3.187930645 0
+id890314 2.099691952 0
+id102874 9.133488079 1
+id200949 7.525406804 0
+id336491 6.798229522 0
+id988766 6.149545358 0
+id21999 1.013546103 1
+id433893 1.282853098 0
+id688932 8.340206657 0
+id394203 3.392345681 1
+id995678 9.587748608 0
+id694339 1.25827049 0
+id256455 5.414840139 0
+id14836 6.177914833 1
+id817128 2.940033056 1
+id803325 9.619267566 0
+id521287 0.320487323 1
+id701472 7.101945547 0
+id850010 6.780876558 1
+id268483 7.903947419 0
+id738781 1.086018804 0
+id28411 1.254350142 1
+id541635 9.057035035 1
+id751101 1.303018406 0
+id826300 6.735484141 1
+id884387 9.097303962 0
+id492414 4.756845946 0
+id268871 4.018189958 0
+id627354 8.874035732 0
+id503932 0.091310535 0
+id163442 3.240808272 0
+id317797 9.050109718 0
+id687857 7.992764162 0
+id871570 2.353719529 0
+id724067 1.617985994 1
+id874076 3.099430273 0
+id927863 3.64738684 0
+id369805 0.973580503 0
+id668376 0.567035975 0
+id717362 6.77104712 1
+id665504 8.956792361 0
+id336637 7.196501733 0
+id60633 5.086557706 0
+id848600 8.270336349 1
+id169514 9.421256192 0
+id690732 9.952887485 0
+id684760 3.361963101 0
+id553502 5.557606017 1
+id214917 6.535326783 0
+id849169 7.903922417 1
+id941921 2.305853849 0
+id784646 3.120883065 1
+id520954 8.756848347 1
+id996355 7.326814334 0
+id96730 6.664821662 0
+id673442 3.081616083 0
+id68305 3.525197877 1
+id653025 7.231586593 0
+id208543 4.622250109 0
+id335725 5.924295793 1
+id980400 2.35961827 0
+id869939 0.736202662 0
+id297563 6.641749611 0
+id852663 5.124175591 0
+id162070 7.910906314 0
+id272875 9.825455416 1
+id163787 2.471002157 0
+id422204 2.915218179 0
+id120197 3.367355865 0
+id33660 4.407010904 1
+id803855 8.755961077 0
+id255048 7.749378138 0
+id690936 7.543593942 0
+id126807 0.251036587 0
+id99016 6.491002725 0
+id883847 1.088929613 0
+id354523 9.618060855 1
+id737255 4.133744506 1
+id990941 5.646270286 0
+id25464 8.272429793 0
+id918375 0.723567779 0
+id537828 4.0147246 0
+id682778 9.130628394 0
+id587547 7.44191196 0
+id670874 4.512303089 0
+id444459 8.723388559 0
+id777456 4.115325335 1
+id452384 2.067203142 0
+id826975 8.348432332 0
+id519567 6.785076383 0
+id84292 7.634117891 0
+id124432 5.099839207 1
+id800145 8.551368147 0
+id153857 3.556301609 0
+id587157 8.758929135 1
+id506262 6.238330685 1
+id634462 0.905623161 0
+id687592 0.290071515 0
+id955526 1.838730964 0
+id181850 0.144560173 1
+id159506 6.045008786 1
+id609051 6.873427527 1
+id963886 0.853758289 0
+id405792 3.909860553 1
+id494172 7.480072965 0
+id964637 3.816858446 0
+id799355 1.226117869 1
+id157111 4.166491482 1
+id114524 8.645968152 0
+id954931 0.158490245 0
+id827034 0.702377072 0
+id689645 2.0658374 0
+id281585 4.640052191 1
+id885624 3.615978697 0
+id577871 2.230068497 0
+id238796 5.161400753 0
+id481035 7.693009686 0
+id972713 6.34411768 1
+id905484 4.55801149 0
+id713511 6.721560405 1
+id512328 4.431789868 0
+id703534 8.974882244 1
+id409904 4.845419164 0
+id577169 1.384047668 0
+id813971 6.785385478 0
+id558483 0.27946972 0
+id892784 1.505769901 0
+id611178 3.97305697 1
+id192732 7.859877791 0
+id917280 1.39759155 0
+id435876 4.928142177 0
+id980722 0.72730016 1
+id308273 3.466919258 1
+id476685 8.925909139 0
+id315883 7.79623385 0
+id935945 3.491058239 1
+id991781 6.268341084 0
+id65199 2.65595714 0
+id226233 6.409789501 0
+id860183 7.567894972 0
+id295209 0.983143735 1
+id544964 1.927387219 0
+id648663 7.914653713 0
+id710165 9.981244413 0
+id392593 7.418579368 0
+id129945 9.344802496 0
+id382621 1.828452452 0
+id901440 9.496967232 0
+id39847 8.962161053 0
+id526460 1.855905599 0
+id477473 1.89081994 1
+id448194 4.814151059 0
+id904184 8.028988991 1
+id747852 3.836260637 1
+id711012 0.341281233 0
+id683879 9.297243459 0
+id789575 5.098990324 0
+id650729 7.941450537 0
+id934302 9.134082902 1
+id555013 8.812444844 1
+id82779 7.657743527 0
+id771444 6.664855358 1
+id821562 5.19036265 0
+id292809 9.909691181 1
+id645690 6.515405399 0
+id223901 6.355198262 1
+id41320 0.397540162 1
+id96181 3.155308448 1
+id147900 4.719584771 0
+id702917 7.72554794 0
+id150640 6.511110099 1
+id518391 9.626758766 0
+id879076 7.61506094 0
+id952031 4.87888217 1
+id10055 5.724749579 0
+id727213 5.731963291 1
+id41961 6.818411434 0
+id257209 9.897638239 1
+id995361 3.151184642 0
+id957918 5.30099154 0
+id975370 9.047264812 0
+id889896 9.578142144 1
+id978164 4.726348119 0
+id90359 6.237920077 0
+id307158 5.49680528 0
+id755940 1.683472511 0
+id995582 1.486566562 0
+id363965 7.723182293 1
+id729124 6.861379688 0
+id871963 1.399526332 1
+id475172 2.745065294 0
+id804699 8.378029915 1
+id625843 4.925798265 1
+id595713 3.273883772 0
+id462604 5.541305429 0
+id106141 9.529928524 0
+id689349 9.243194854 0
+id639003 6.594759182 1
+id393896 8.69916317 1
+id450307 7.291170542 1
Added: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh (rev 0)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh 2015-03-17 13:54:45 UTC (rev 1931)
@@ -0,0 +1,82 @@
+#!/bin/bash
+# This script runs checks on ProbABEL's pacoxph module
+# These are almost identical to the checks in test_cox.sh, but this
+# time without covariate data (see bug #1266).
+
+
+echo "Analysing Cox model without covariates..."
+
+scriptdir=$(dirname $0)
+
+if [ -z ${PA_BINDIR} ]; then
+ PA_BINDIR="${scriptdir}/../src/"
+fi
+if [ -z ${srcdir} ]; then
+ srcdir="."
+ PA_BINDIR=${scriptdir}/../src/
+fi
+
+. ${scriptdir}/run_diff.sh
+
+inputdir=${scriptdir}/inputfiles
+pacoxph=${PA_BINDIR}/pacoxph
+
+# Redirect all output to file descriptor 3 to /dev/null except if
+# the first argument is "verbose" then redirect handle 3 to stdout
+exec 3>/dev/null
+if [ "$1" = "verbose" ]; then
+ echo "Verbose mode ON"
+ exec 3>&1
+fi
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.mldose \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ -c 19 \
+ -o coxph_dose_nocovar \
+ >& 3
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.dose.fvi \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ -c 19 \
+ -o coxph_dose_nocovar_fv \
+ >& 3
+
+run_diff coxph_dose_nocovar_add.out.txt coxph_dose_nocovar_fv_add.out.txt \
+ "pacoxph check: dose vs. dose_fv"
+
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.mlprob \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob_nocovar \
+ >& 3
+
+run_diff coxph_dose_nocovar_add.out.txt coxph_prob_nocovar_add.out.txt \
+ "pacoxph check: dose vs. prob" -I SNP
+
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.prob.fvi \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob_nocovar_fv \
+ >& 3
+
+for model in add domin recess over_domin 2df; do
+ run_diff coxph_prob_nocovar_${model}.out.txt \
+ coxph_prob_nocovar_fv_${model}.out.txt \
+ "pacoxph check ($model model): prob vs. prob_fv"
+done
Property changes on: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh
___________________________________________________________________
Added: svn:executable
+ *
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/coxph_data.cpp 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/coxph_data.cpp 2015-03-17 13:54:45 UTC (rev 1931)
@@ -5,7 +5,7 @@
* Author: mkooyman
*
*
- * Copyright (C) 2009--2014 Various members of the GenABEL team. See
+ * Copyright (C) 2009--2015 Various members of the GenABEL team. See
* the SVN commit logs for more details.
*
* This program is free software; you can redistribute it and/or
@@ -116,16 +116,16 @@
exit(1);
}
- X.reinit(nids, ncov);
+ X.reinit(nids, (ncov + 1));
stime.reinit(nids, 1);
sstat.reinit(nids, 1);
weights.reinit(nids, 1);
offset.reinit(nids, 1);
strata.reinit(nids, 1);
order.reinit(nids, 1);
+
for (int i = 0; i < nids; i++)
{
- // X.put(1.,i,0);
stime[i] = (phed.Y).get(i, 0);
sstat[i] = static_cast<int>((phed.Y).get(i, 1));
if (sstat[i] != 1 && sstat[i] != 0)
@@ -137,12 +137,21 @@
}
}
- for (int j = 0; j < phed.ncov; j++)
+ // Add a column with a constant (=1) to the X matrix (the mean)
+ for (int i = 0; i < nids; i++)
{
+ X.put(1., i, 0);
+ }
+
+ for (int j = 1; j <= phed.ncov; j++)
+ {
for (int i = 0; i < nids; i++)
- X.put((phed.X).get(i, j), i, j);
+ {
+ X.put((phed.X).get(i, j - 1), i, j);
+ }
}
+
if (snpnum > 0)
{
for (int j = 0; j < ngpreds; j++)
@@ -246,7 +255,7 @@
gend->get_var(snpnum * ngpreds + j, snpdata);
for (int i = 0; i < nids; i++) {
- X.put(snpdata[i], (ncov - j - 1), order[i]);
+ X.put(snpdata[i], (ncov - j), order[i]);
if (std::isnan(snpdata[i])) {
masked_data[order[i]] = 1;
// snp not masked
@@ -398,8 +407,6 @@
(cdata.offset).column_mean(0);
mematrix<double> means(X.nrow, 1);
- int maxiterinput = maxiter;
-
for (int i = 0; i < X.nrow; i++)
{
beta[i] = 0.;
@@ -418,6 +425,10 @@
// R's coxph()
double sctest = 1.0;
+ // Set the maximum number of iterations that coxfit2() will run to
+ // the default value from the class definition.
+ int maxiterinput = maxiter;
+
// When using Eigen coxfit2 needs to be called in a slightly
// different way (i.e. the .data()-part needs to be added).
#if EIGEN
@@ -449,18 +460,39 @@
// to NA if flag < nvar (with nvar = ncol(x)) and maxiter >
// 0. coxph.R then checks for any NAs in the coefficients and
// outputs the warning message if NAs were found.
- if (flag < X.nrow && maxiter > 0) {
- cerr << "Warning for " << snpinfo.name[cursnp]
- << ": X matrix deemed to be singular,"
- << " setting beta and se to 'NaN'\n";
- setToNAN = true;
+ if (flag < X.nrow)
+ {
+#if EIGEN
+ int which_sing = NAN;
+ MatrixXd imateigen = imat.data;
+ VectorXd imatdiag = imateigen.diagonal();
+
+ // Start at i=1 to exclude the beta coefficient for the
+ // (constant) mean from the check.
+ for (int i=1; i < imatdiag.size(); i++)
+ {
+ if (imatdiag[i] == 0)
+ {
+ which_sing = i;
+ setToNAN = true;
+ std::cerr << "Warning for " << snpinfo.name[cursnp]
+ << ": X matrix deemed to be singular (variable "
+ << which_sing + 1 << ")" << std::endl;
+ }
+ }
+#else
+ std::cerr << "Warning for " << snpinfo.name[cursnp]
+ << ": can't check for singular X matrix."
+ << " Please compile ProbABEL with Eigen support to fix this."
+ << std::endl;
+#endif
}
if (niter >= maxiterinput)
{
cerr << "Warning for " << snpinfo.name[cursnp]
- << ": nr of iterations > MAXITER (" << maxiterinput << "): "
+ << ": nr of iterations >= MAXITER (" << maxiterinput << "): "
<< niter << endl;
}
@@ -485,7 +517,8 @@
// (incl. covariates), maybe stick to only checking the SNP
// coefficient?
for (int i = 0; i < infs.size(); i++) {
- if (infs[i] > eps && infs[i] > sqrt(eps) * abs(betaeigen[i])) {
+ if (infs[i] > eps &&
+ infs[i] > sqrt(eps) * abs(betaeigen[i])) {
problems = true;
}
}
@@ -497,9 +530,6 @@
setToNAN = true;
}
-
- // cout << "beta values for SNP: " << snpinfo.name[cursnp]
- // << " are: " << betaeigen << std::endl;
#else
cerr << "Warning for " << snpinfo.name[cursnp]
<< ": can't check for infinite betas."
Modified: branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/phedata.cpp
===================================================================
--- branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/phedata.cpp 2015-03-16 22:01:14 UTC (rev 1930)
+++ branches/ProbABEL-v0.4.4-hotfix/ProbABEL/src/phedata.cpp 2015-03-17 13:54:45 UTC (rev 1931)
@@ -135,17 +135,14 @@
model = model + tmp;
}
n_model_terms = 0;
-#if COXPH
- model = model + " ) ~ ";
-#else
- model = model + " ) ~ mu + ";
+
+ model = model + " ) ~ mu";
model_terms[n_model_terms++] = "mu";
-#endif
if (nphenocols > noutcomes + 1)
{
infile >> tmp;
- model = model + tmp;
+ model = model + " + " + tmp;
model_terms[n_model_terms++] = tmp;
for (int i = (2 + noutcomes); i < nphenocols; i++)
{
More information about the Genabel-commits
mailing list