snp.logistic(data, response.var, snp.var, main.vars=NULL, int.vars=NULL, strata.var=NULL, op=NULL)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). genetic.model, reltol, maxiter, and
optimizer (see details). The default is NULL.
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\frac{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|S=s) = (1 - f_s)^2$$ $$P(G = 1|S=s) = 2f_s(1 - f_s)$$ $$P(G = 2|S=s) = f_s^2.$$ If $xi_s = log(f_S/(1 - f_s))$, then $$\log \left( \frac{P(G = 1)}{P(G = 0)} \right) = \log(2) + \xi_s$$ and $$\log \left( \frac{P(G = 2)}{P(G = 0)} \right) = 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) = \frac{\exp(\theta_{s}(d, g|Z))}{\sum_{d,g} \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 \left( \frac{\phi^{2}(\phi^{2} - \psi^2)}{(\phi^{2} + \psi^2)^2} \right)$$ 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".
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))
Run the code above in your browser using DataLab