POM-PS (proportional odds model adjusted for propensity score) is a novel statistical approach for testing association of one or more secondary traits with a single genetic marker within the framework of a case-control design. The R function pomps implements this association test. For details of this statistical method, please refer/cite:

Ray, D. and Basu, S. "A Novel Association Test for Multiple Secondary Phenotypes from a Case-Control GWAS". Genetic Epidemiology, 41:413-426, 2017. PMID: 28393390

Key Words: Case-control design; Cross-phenotype association; GWAS; Joint modeling; Multiple traits; Multivariate analysis; Propensity Score; Proportional Odds Model; Secondary traits; Stratification Score

R (>= 2.14.0), MASS, survey

require(devtools)
source_url("https://github.com/RayDebashree/POM-PS/blob/master/pomps_v1.15.R?raw=TRUE")

It is recommended to download/copy the stand-alone R program in this repository, save it in your local directory of choice and source() it from your local directory. When a new version of the software is available, older versions may be removed from this repository, and the above devtools::source_url() technique may not work.

Version 1.15 - February 28, 2017

First public release of the software.

pomps(Y, X, D=NULL, COV=NULL, PS=NULL, method="POM-PS", add.D.as.COV=FALSE, msg.mute=FALSE, no.format.check=FALSE, ...)

source("pomps_v1.15.R")
# simulate 2 phenotypes on 1000 individuals
Y<-mvrnorm(n=1000, mu=c(0,0), Sigma=matrix(c(1,0.2,0.2,1),2,2))

# simulate 3 covariates for disease status
Sig<-matrix(c(1,0.1,-0.45,0.1,1,0.6,-0.45,0.6,1),3,3)
Z<-mvrnorm(n=1000, mu=c(0,0,0), Sigma=Sig)

# simulate error for disease status
E<-rt(n=1000, df=1,ncp=0)

# simulate disease status
D<-rep(0,1000)
D[which(rowSums(cbind(Y,Z,E))>=0)]<-1

# simulate a single marker for 1000 individuals
X<-matrix(rbinom(n=1000, size=2, prob=0.2), ncol=1)	# additive model

# required data-frame formats
Y<-as.data.frame(Y)
D<-as.data.frame(D)
Z<-as.data.frame(Z)
X<-as.data.frame(X)

# unique column names for the data-frames
colnames(Y)<-paste("Y",1:2,sep="")
colnames(Z)<-paste("Z",1:3,sep="")
colnames(X)<-"X"
colnames(D)<-"D"

## apply POM to test association
out1<-pomps(Y=Y, X=X, COV=NULL, method="POM")
# POM test statistic and p-value
t1<-out1$stat
p1<-out1$pvalue

## apply POM-PS to test association
S=getPS(Y=Y, D=D, covars=Z)
out2<-pomps(Y=Y, X=X, COV=NULL, PS=S)
# POM-PS test statistic and p-value
t2<-out2$stat
p2<-out2$pvalue

1. The proportional odds framework of POM-PS allows one or more secondary phenotypes, which may be binary and/or continuous.

2. The method POM-PS and its software is designed for unrelated individuals. If you have two cohorts with overlapping samples and you want to analyze the combined sample, it is desirable to exclude the overlapping individuals, and any related individuals.

3. The genotype X takes values 0, 1 or 2. For a given sample, if X has only two possible values, a logistic regression (glm) is used instead of proportional odds regression (polr).

4. POM-PS is a single-variant association test, and is not expected to work well for rare variants (i.e., when the variant allele-frequency is very low).

5. Apart from function pomps(), this software has other relevant functions such as getPS() (for calculating propensity score of an individual given secondary phenotypes, case-control status and other covariates - excluding genetic variants) and pom.ipw() (for implementing an inverse-probability-weighted proportional odds regression model). Details of these functions are provided in the attached manual.

   • Remember: The estimated propensity score using the getPS() function needs to be obtained only once for a given GWAS.
   • Caution: While pomps() removes missing observations (if any) from the input data, getPS() is not programmed to handle missing observations (NA values) and hence we need to provide complete (non-missing) data as input for getPS().
   • Caution: pom.ipw() often encounters convergence issues; use with caution.

6. POM-PS currently analyzes secondary phenotypes when the primary trait is binary. For primary trait with three or more ordinal categories, one may employ a proportional odds model for calculating the propensity score.