PMR_summary {PPMR}

PMR model with summary data

Description

Fit the probabilistic MR model with summary data while accounting for the correlated instruments and horizontal pleiotropy in TWAS framework.

Usage

PMR_summary(betaxin, betayin, Sigma1sin, Sigma2sin, samplen1, samplen2,
  gammain, alphain, max_iterin, epsin)

Arguments

betaxin

the cis-SNP effect size vector for one specific gene in eQTL data, which must be calculated based on both the standardized gene expression value and the standardized cis-genotype matrix.

betayin

the cis-SNP effect size vector for one specific gene in GWAS data, which be calculated based on both the standardized complex trait value and the standardized cis-genotype matrix.

Sigma1sin

the LD matrix in eQTL data.

Sigma2sin

the LD matrix in GWAS data. Both Sigma2sin and sigma1sin are often the same from the reference panel.

samplen1

the sample size of eQTL data.

samplen2

the sample size of GWAS data.

gammain

indicator variable for constrained model, with 1 for the null model that there is no horizontal pleiotropy.

alphain

indicator variable for constrained model, with 1 for the null model that there is no causal effect.

max_iterin

The maximum iteration.

epsin

The convergence tolerance of the absolute value of the difference between the nth and (n+1)th log likelihood.

Value

A list of estimates of model parameters, including the causal effect alpha, the horizontal pleiotropy effect gamma, and the two corresponding p values.

Author(s)

Zhongshang Yuan, Xiang Zhou.

Examples

data(Examplesummary)
attach(Examplesummary)
fmH1=PMR_summary(betaxin=betax,betayin=betay,Sigma1sin=Sigma1,Sigma2sin=Sigma2,
      samplen1=n1,samplen2=n2,gammain=0,alphain=0,max_iterin =1000, epsin=1e-5)
fmH0alpha=PMR_summary(betaxin=betax,betayin=betay,Sigma1sin=Sigma1,Sigma2sin=Sigma2,
           samplen1=n1,samplen2=n2,gammain=0,alphain=1,max_iterin =1000, epsin=1e-5)
fmH0gamma=PMR_summary(betaxin=betax,betayin=betay,Sigma1sin=Sigma1,Sigma2sin=Sigma2,
           samplen1=n1,samplen2=n2,gammain=1,alphain=0,max_iterin =1000, epsin=1e-5)
loglikH1=max(fmH1$loglik,na.rm=TRUE)
loglikH0alpha=max(fmH0alpha$loglik,na.rm=TRUE)
loglikH0gamma=max(fmH0gamma$loglik,na.rm=TRUE)
stat_alpha = 2 * (loglikH1 - loglikH0alpha)
pvalue_alpha = pchisq(stat_alpha,1,lower.tail=FALSE)
stat_gamma = 2 * (loglikH1 - loglikH0gamma)
pvalue_gamma = pchisq(stat_gamma,1,lower.tail=FALSE)