Interaction G*E effect/SE/ LRT test statistics incorrect

Questions about ProbABEL are welcome here.
Forum rules
Please remember not to post any sensitive data on this public forum.
The first few posts of newly registered users will be moderated in order to filter out any spammers.

When get a solution to the problem you posted, please change the topic name (e.g. from "how to ..." to "[SOLVED] how to ..."). This will make it easier for the community to follow the posts yet to be attended.
ericyin
Posts: 4
Joined: Mon Jul 20, 2015 10:03 am

Interaction G*E effect/SE/ LRT test statistics incorrect

Postby ericyin » Tue Sep 22, 2015 5:07 pm

Hi,

Two problems I found for interaction effect in both palogist and pacoxph.

I was firstly running a logistic regression with interaction effect (drug * SNP). I found the LRT test statistics is different with results using R. (The estimate and standard error is the same!)

I am running coxph with interaction effect (also drug *SNP). I found the pacoxph output for interaction is incorrect (estimate, se, LRT test statistics) with results using R (survival).

See outputs here: (I am running Statins*SNP interaction here)

R outputs

Code: Select all

coxph(formula = Surv(Days_FU1, primary_endpoint) ~ as.numeric(age) +
    as.numeric(BMI) + CRF + PCI + PAST_MI + CABG + ACEI_pre +
    Beta_blocker_pre + Statins + Beta_blocker + Aldosterone +
    Diabetes + snp3_dose + Statins * snp3_dose, data = dataset)

  n= 1274, number of events= 182
   (83 observations deleted due to missingness)

                       coef exp(coef)  se(coef)      z Pr(>|z|)
as.numeric(age)    0.041705  1.042587  0.007617  5.475 4.37e-08 ***
as.numeric(BMI)    0.015886  1.016013  0.010544  1.507  0.13193
CRF                0.355794  1.427313  0.208758  1.704  0.08832 .
PCI               -0.482626  0.617161  0.197648 -2.442  0.01461 *
PAST_MI            0.520320  1.682566  0.163494  3.183  0.00146 **
CABG              -0.252058  0.777200  0.138693 -1.817  0.06916 .
ACEI_pre           0.319004  1.375757  0.158544  2.012  0.04421 *
Beta_blocker_pre   0.409746  1.506434  0.170574  2.402  0.01630 *
Statins           -0.924099  0.396889  0.976754 -0.946  0.34410
Beta_blocker      -0.501741  0.605476  0.187847 -2.671  0.00756 **
Aldosterone        0.521484  1.684526  0.254144  2.052  0.04018 *
Diabetes           0.239759  1.270943  0.168461  1.423  0.15467
snp3_dose          0.444876  1.560297  0.512900  0.867  0.38574
Statins:snp3_dose  0.260580  1.297682  0.535972  0.486  0.62684


pacoxph output:

name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele beta_mu sebeta_mu beta_age sebeta_age beta_BMI sebeta_BMI beta_CRF sebeta_CRF beta_PCI sebeta_PCI beta_CABG sebeta_CABG beta_PAST_MI sebeta_PAST_MI beta_Beta_blocker_pre sebeta_Beta_blocker_pre beta_ACEI_pre sebeta_ACEI_pre beta_Diabetes sebeta_Diabetes beta_Beta_blocker sebeta_Beta_blocker beta_Statins sebeta_Statins beta_Aldosterone sebeta_Aldosterone beta_SNP_add sebeta_SNP_add beta_SNP_Statins sebeta_SNP_Statins chi2_SNP

rs300251 C T 0.75383 0.24617 0.99596 0.1 1274 0.754266 0 0 0.0415631 0.00762207 0.016612 0.010929 0.353853 0.208645 -0.483521 0.197588 -0.254933 0.138657 0.513987 0.163178 0.409519 0.170432 0.313767 0.158361 0.240251 0.168653 -0.499339 0.187663 0.516968 0.254194 0.896249 0.196951 -0.223234 0.142819 23.0213



As we see here, the last five columns of outputs

beta_SNP_add sebeta_SNP_add beta_SNP_Statins sebeta_SNP_Statins chi2_SNP

are different from R output...




Also there is a small bug in the --interaction in pacoxph. See below:

Code: Select all

pacoxph -p ProbABEL_pheno_coxph_Statins_1.txt -d ACSchr_5.dose.fvi -i 5_ACS_QC.gen.Status.out.mlinfo --allcov --interaction=12 --out ACS_DC_5_coxph_Statins_int_alloutput


Check the file
ProbABEL_pheno_coxph_Statins_1.txt
ID_2 Days_FU1 primary_endpoint age BMI CRF PCI CABG PAST_MI Beta_blocker_pre ACEI_pre Diabetes Beta_blocker Statins Aldosterone


“Statins” is the 11th covariate, but you need to state --interaction=12 to run Statins*SNP (I think the programmer may forget Survival analysis (pacoxph) using two columns to save phenotype data while logistic (palogist) using 1 column.)

Code: Select all

probabel v. 0.4.5
(C) Yurii Aulchenko, Lennart C. Karssen, Maarten Kooyman, Maksim Struchalin, The
 GenABEL team, EMC Rotterdam

Using EIGEN version 3.2.0 for matrix operations
Options in effect:
         --pheno       = ProbABEL_pheno_coxph_Statins_1.txt
         --info        = 5_ACS_QC.gen.Status.out.mlinfo
         --dose        = ACSchr_5.dose.fvi
                     (using FVF data)
         --ntraits     = 2
         --ngpreds     = 1
         --interaction = 12
         --interaction_only = 0
         --mmscore = not in output
         --map     = not in output
         --nids    = estimated from data
         --chrom   = not in output
         --out     = ACS_DC_5_coxph_Statins_int_alloutput
         --skipd   = 2
         --separat = " "
         --score   = OFF
         --nohead  = OFF
         --allcov  = ON
         --robust  = OFF

Reading info data...
Number of SNPs = 1570641
Reading phenotype data...
Actual number of people in phenofile = 1357; using all of these
Coxph model: ( Days_FU1 , primary_endpoint ) ~ mu + age + BMI + CRF + PCI + CABG
 + PAST_MI + Beta_blocker_pre + ACEI_pre + Diabetes + Beta_blocker + Statins + A
ldosterone + SNP_A1 + Statins*SNP_A1
Reading genotype data... done.
 loaded null data... estimated null model... formed regression object...
Analysis:  0.19%...

Return to “ProbABEL”

Who is online

Users browsing this forum: No registered users and 1 guest