search

suchestoncampbelllab / gwasurvivr

GWAS Survival Package in R
11 stars 12 forks source link

gwasurvivr and R "survival" package give different results #7

Closed Qiaolan closed 4 years ago

Qiaolan commented 4 years ago

Hi,

I find that given the same genotype data and phenotype data. gwasurvivr and R survival package cox.ph() give different p-value (both significant, but p-values are quite different, for example, one might be 1e-6 and the other 1e-7)and coefs.

Could you please give me any help? Is there any possible reason?

Thanks!

aarizvi commented 4 years ago

Please send a reproducible example and I can look at it with you. In theory they shouldn't be off by an order of magnitude unless the genotyping (or imputation) quality is poor. We are using the same survival:::coxph.fit internal function in this package so the estimators are the same. The only difference is that in gwasurvivr we put in initial estimates that we obtain by fitting the model just the covariates first.

Qiaolan commented 4 years ago

Thanks for your reply. I am sorry the data is in the server I cannot copy it.

Below are codes of two methods

  1. R survival::coxph: model <- coxph(Surv(age_onset,status) ~ cov_match[,28+i] +sex+pc1+pc2+pc3+pc4+pc5,data=cov_match,ties = "efron")

  2. gwasurvivr:

bed.file <- "data_clean/Two/cov_id.bed" cov <- read.table("data_clean/Two/cov_match.txt",header=T) sample.ids <- as.character(cov$merged_id)

plinkCoxSurv(bed.file=bed.file, covariate.file=cov, id.column="merged_id", sample.ids=sample.ids, time.to.event="age_onset", event="status", covariates=c("sex", "pc1","pc2","pc3","pc4","pc5"), inter.term=NULL, print.covs="only", out.file= "/fs/project/PAS1501/gaolab/qldeng/Eden/data_clean/Two/result", chunk.size=100, maf.filter=0.01, flip.dosage=TRUE, verbose=TRUE, clusterObj=cl)

I transferred bed file to 0/1/2 and it is contained in cov_match.txt used by cox.ph()

Also, I checked the source code of gwasurvivr and I did find that it uses coxph.fit() rather than coxph. Is there any difference? I find few manual about coxph.fit()

At last, you mentioned "The only difference is that in gwasurvivr we put in initial estimates that we obtain by fitting the model just the covariates first." When I use coxph() in R survial, is there anything I need to modify?

Thanks for your help!

Qiaolan commented 4 years ago

Hi and thanks for your time!

Today I checked the sample data within the gwasurvivr package. The results of survival::coxph and gwasurvivr match. Moreover, I find there's no difference if I set initial values for coxph.fit, because both coxph and coxph.fit with initials gave me the same coefs. Is this possible?

At last, returning to my original problem, my data is case-only, that is, all subjects experienced the event. Is this a possible reason that two methods gave me different coefs?

Thank you!