snp.logistic: Logistic regression analysis for a single SNP
In CGEN: An R package for analysis of case-control studies in genetic epidemiology

Description Usage Arguments Details Value References See Also Examples

Performs logistic regression including a particular SNP (G) and a set of covariates (X) that could include environmental covariates or/and other genetic variables. Included are three analysis options: (i) Unconstrained maximum-likelihood: This method is equivalent to prospective logistic regression analysis and corresponds to maximum-likelihood analysis of case-control data allowing the joint distribution of all the factors (the SNP of interest and all other covariates) of the model to be completely unrestricted (non-parametric) (ii) Constrained maximum-likelihood: This method performs maximum-likelihood analysis of case-control data under the assumption of HWE and indepenence between the SNP and other factors of the model.The analysis allows the assumptions of HWE and independence to be valid only conditional on certain stratification variables (S), such as self reported ethnicity or principal compoenets of population stratification. (iii) Empirical-Bayes: This method uses an empirical-Bayes type "shrinkage estimation" technique to trade-off bias and variance between the constrained and unconstrained maximum-likelihood estimators.

1 2	snp.logistic(data, response.var, snp.var, main.vars=NULL, int.vars=NULL, strata.var=NULL, op=NULL)

`data`	Data frame containing all the data. No default.
`response.var`	Name of the binary response variable coded as 0 (controls) and 1 (cases). No default.
`snp.var`	Name of the SNP variable, which must be coded 0-1-2 (or 0-1). The SNP will be included as a main effect in the model. No default.
`main.vars`	Character vector of variable names or a formula for all covariates of interest which need to be included in the model as main effects. The default is NULL, so that only the SNP variable will be included as a main effect in the model.
`int.vars`	Character vector of variable names or a formula for all covariates of interest that will interact with the SNP variable. The default is NULL, so that no interactions will be in the model.
`strata.var`	Name of the stratification variable or a formula (see details for more info). If `strata.var`="SVAR", where "SVAR" is a factor or character variable in `data`, then "SVAR" will be treated as categorical. Otherwise, "SVAR" is treated as a continuous variable. The default is NULL (1 stratum).

`op`	A list with names `genetic.model`, `reltol`, `maxiter`, and `optimizer` (see details). The default is NULL.

Note: Non-dummy continuous variables should be scaled for stability of the algorithm. The scale function can be used for this.

The data is first fit using standard logistic regression. The estimated parameters from the standard logistic regression are then used as the initial estimates for the constrained model. For this, the optim() function is used to compute the maximum likelihood estimates and the estimated covariance matrix. The empirical Bayes estimates are then computed by combining both sets of estimated parameters (see below). The "strata" option, that is relevent for the CML and EB method, allows the assumption of HWE and G-X independence to be valid only conditional on a given set of other factors. If a single categorical variable name is provided, then the unique levels of the variable will be used to define categorical strata. Otherwise it is assumed that strata.var defines a parametric model for variation of allele frequency of the SNP as a function of the variables included. No assumption is made about the relationship between X and S. Typically, S would include self reported ethnicity, study, center/geographic region and principal components of population stratification. The CML method with the "strata" defined by principal compoenents of population stratification can be viewed as a generalization of adjusted case-only method described in Bhattacharjee et al. (2010). More details of the individual methods follow.

Definition of the likelihood under the gene-environment independence assumption:

Let D = 0, 1 be the case-control status, G = 0, 1, 2 denote the SNP genotype, S denote the stratification variable(s) and X denote the set of all other factors to be included in the regression model. Suppose the risk of the disease (D), given G, X and S can be described by a logistic regression model of the form

log(Pr(D=1)/Pr(D=0)) = alpha + Z*beta

where Z is the entire design matrix (including G, X, possibly S and their interaction with X) and beta is the vector of associated regression coefficients. The CML method assumes Pr(G|X,S)=Pr(G|S), i.e., G and X are conditionally independent given S. The current implementation of the CML method also assume the SNP genotype frequency follows HWE given S=s, although this is not necessary in general. Thus, if f_s denotes the allele frequency given S=s, then

P(G = 0) = (1 - f_s)^2,

P(G = 1) = 2*f_s*(1 - f_s),

P(G = 2) = (f_s)^2.

If xi_s = log(f_S/(1 - f_s)), then

log(P(G=1)/P(G=0)) = log(2) + xi_s

and

log(P(G=2)/P(G=0)) = 2*xi_s.

Chatterjee and Carroll (2005) showed that under the above constraints, the maximum-likelihood estimate for the beta coefficients under case-control design can be obtained based on a simple conditional likelihood of the form

P`(D=d, G=g | Z, S) = exp(theta_s(d,g|Z))/SUM[exp(theta_s(d,g|Z))]

where the sum is taken over the 6 combinations of d and g and theta_s(d,g) = d*alpha` + d*Z*beta + I(g=1)*log(2) + g*xi_s. If S is a single categorical variable, then a separate xi_s is allowed for each S=s. Otherwise it is assumed xi_s=V_s*gamma, where V_s is the design matrix associated with the stratification and gamma is the vector of stratification parameters. If for example, S is specified as "strata=~PC1+PC2+...PCK" where PCk's denote principal components of population stratification, then it is assumed that the allele frequency of the SNP varies in directions of the different principal components in a logistic linear fashion.

Definition of the empirical bayes estimates:
Let beta_UML be the parameter estimates from standard logistic regression, and let eta = (beta_CML, xi_CML) be the estimates under the gene-environment independence assumption. Let psi = beta_UML - beta_CML, and phi^2 be the vector of variances of beta_UML. Define diagonal matrices of weights to be W1 = diag(psi^2/(psi^2 + phi^2)) and W2 = diag(phi^2/(psi^2 + phi^2)), where psi^2 is the elementwise product of the vector psi. Now, the empirical bayes parameter estimates are

beta_EB = W1*beta_UML + W2*beta_CML.

For the estimated covariance matrix, define the diagonal matrix

A = diag(phi^2*(phi^2 - psi^2)/(phi^2 + psi^2)^2),

where again the exponentiation is the elementwise product of the vectors. If I is the pxp identity matrix and we define the px2p matrix C = (A, I - A), then the estimated covariance matrix is

VAR(beta_EB) = C*COV(beta_UML, beta_CML)*C'.

The covariance term COV(beta_UML, beta_CML) is obtained using an influence function method (see Chen YH, Chatterjee N, and Carroll R. for details about the above formulation of the empirical-Bayes method).

Options list:

Below are the names for the options list op. All names have default values if they are not specified.

genetic.model 0-3: The genetic model for the SNP. 0=additive, 1=dominant, 2=recessive, 3=general (co-dominant).
reltol Stopping tolerance. The default is 1e-8.
maxiter Maximum number of iterations. The default is 100.
optimizer One of "BFGS", "CG", "L-BFGS-B", "Nelder-Mead", "SANN". The default is "BFGS".

A list containing sublists with names UML (unconstrained maximum likelihood), CML (constrained maximum likelihood), and EB (empirical Bayes). Each sublist contains the parameter estimates (parms) and covariance matrix (cov). The lists UML and CML also contain the log-likelihood (loglike). The list CML also contains the results for the stratum specific allele frequencies under the HWE assumption (strata.parms and strata.cov). The EB sublist contains the joint UML-CML covariance matrix.

Mukherjee B, Chatterjee N. Exploiting gene-environment independence in analysis of case-control studies: An empirical Bayes approach to trade-off between bias and efficiency. Biometrics 2008, 64(3):685-94.

Mukherjee B et al. Tests for gene-environment interaction from case-control data: a novel study of type I error, power and designs. Genetic Epidemiology, 2008, 32:615-26.

Chatterjee, N. and Carroll, R. Semiparametric maximum likelihood estimation exploting gene-environment independence in case-control studies. Biometrika, 2005, 92, 2, pp.399-418.

Chen YH, Chatterjee N, Carroll R. Shrinkage estimators for robust and efficient inference in haplotype-based case-control studies. Journal of the American Statistical Association, 2009, 104: 220-233.

Bhattacharjee S, Wang Z, Ciampa J, Kraft P, Chanock S, Yu K, Chatterjee N Using Principal Components of Genetic Variation for Robust and Powerful Detection of Gene-Gene Interactions in Case-Control and Case-Only studies. American Journal of Human Genetics, 2010, 86(3):331-342.

snp.score, snp.matched

 # Use the ovarian cancer data
 data(Xdata, package="CGEN")

 # Fit using a stratification (categorical) variable
 ret <- snp.logistic(Xdata, "case.control", "BRCA.status",
                     main.vars=c("oral.years", "n.children"), 
                     int.vars=c("oral.years", "n.children"), 
                     strata.var=~factor(ethnic.group))

 # Compute a summary table for the models
 getSummary(ret)

 # Compute a Wald test for the main effect of the SNP and interaction
 getWaldTest(ret, c("BRCA.status", "BRCA.status:oral.years", "BRCA.status:n.children"))

 # Fit the same model as above using formulas
 ret2 <- snp.logistic(Xdata, "case.control", "BRCA.status",
                     main.vars=~oral.years + n.children, 
                     int.vars=~oral.years + n.children, 
                     strata.var=~factor(ethnic.group))

Loading required package: survival
Loading required package: mvtnorm
Warning message:
In snp.logistic(Xdata, "case.control", "BRCA.status", main.vars = c("oral.years",  :
  For stability of the algorithm, continuous variables should be scaled
$UML
                          Estimate  Std.Error    Z.value       Pvalue
Intercept              -0.10879020 0.09785000 -1.1118059 2.662216e-01
oral.years             -0.06016844 0.02558151 -2.3520282 1.867136e-02
n.children             -0.02340819 0.02859525 -0.8186044 4.130122e-01
BRCA.status             3.68693982 0.65688833  5.6127346 1.991540e-08
BRCA.status:oral.years  0.05714218 0.10114601  0.5649474 5.721095e-01
BRCA.status:n.children -0.20527943 0.21238273 -0.9665543 3.337669e-01

$CML
                          Estimate  Std.Error    Z.value       Pvalue
Intercept              -0.10406760 0.09852661 -1.0562385 2.908593e-01
oral.years             -0.06023353 0.02539981 -2.3714166 1.772005e-02
n.children             -0.02505173 0.02891007 -0.8665401 3.861940e-01
BRCA.status             3.24287591 0.32650820  9.9319891 3.021666e-23
BRCA.status:oral.years  0.08376709 0.03136942  2.6703423 7.577395e-03
BRCA.status:n.children -0.04344904 0.05139541 -0.8453875 3.978945e-01

$EB
                          Estimate  Std.Error    Z.value       Pvalue
Intercept              -0.10407857 0.09854088 -1.0561969 2.908783e-01
oral.years             -0.06023353 0.02539993 -2.3714053 1.772059e-02
n.children             -0.02504632 0.02892130 -0.8660165 3.864811e-01
BRCA.status             3.38215836 0.53517626  6.3197092 2.620559e-10
BRCA.status:oral.years  0.08204177 0.03668498  2.2363856 2.532652e-02
BRCA.status:n.children -0.10289443 0.17905124 -0.5746647 5.655181e-01

$UML
$UML$test
[1] 112.8177

$UML$df
[1] 3

$UML$pvalue
[1] 2.715665e-24


$CML
$CML$test
[1] 121.0553

$CML$df
[1] 3

$CML$pvalue
[1] 4.572369e-26


$EB
$EB$test
[1] 119.0982

$EB$df
[1] 3

$EB$pvalue
[1] 1.206813e-25


Warning message:
In snp.logistic(Xdata, "case.control", "BRCA.status", main.vars = ~oral.years +  :
  For stability of the algorithm, continuous variables should be scaled