Improving the power of genetic association tests with imperfect phenotype derived from electronic medical records

Abstract

To reduce costs and improve clinical relevance of genetic studies, there has been increasing interest in performing such studies in hospital-based cohorts by linking phenotypes extracted from electronic medical records (EMRs) to genotypes assessed in routinely collected medical samples. A fundamental difficulty in implementing such studies is extracting accurate information about disease outcomes and important clinical covariates from large numbers of EMRs. Recently, numerous algorithms have been developed to infer phenotypes by combining information from multiple structured and unstructured variables extracted from EMRs. Although these algorithms are quite accurate, they typically do not provide perfect classification due to the difficulty in inferring meaning from the text. Some algorithms can produce for each patient a probability that the patient is a disease case. This probability can be thresholded to define case–control status, and this estimated case–control status has been used to replicate known genetic associations in EMR-based studies. However, using the estimated disease status in place of true disease status results in outcome misclassification, which can diminish test power and bias odds ratio estimates. We propose to instead directly model the algorithm-derived probability of being a case. We demonstrate how our approach improves test power and effect estimation in simulation studies, and we describe its performance in a study of rheumatoid arthritis. Our work provides an easily implemented solution to a major practical challenge that arises in the use of EMR data, which can facilitate the use of EMR infrastructure for more powerful, cost-effective, and diverse genetic studies.

Appendix

Design A

In Design A, we consider a random sample of size $n$ from the entire EMR data and calculate $\hat{p}_{D}$ for everyone in the sample. Using assumption (${\fancyscript{A}}$), we see $P(\hat{p}_{D}> c \mid {\mathbf{X}}) = P(\hat{p}_{D}> c \mid D=1) g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}) + P(\hat{p}_{D}> c \mid D=0)(1-g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})).$ Then, since for any positive random variable $T$, $E[T] = \int _0^\infty P(T > c)dc$, we have: $E[\hat{p}_{D}\mid {\mathbf{X}}] = \int _0^1P(\hat{p}_{D}>c \mid {\mathbf{X}})dc = g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}) \int _0^1 P(\hat{p}_{D}> c \mid D=1)dc + (1-g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}))\int _0^1P(\hat{p}_{D}> c \mid D=0)dc = \zeta _0 + (\zeta _1 - \zeta _0)g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})$ where $\zeta _d = E[\hat{p}_{D}\mid D=d].$ Thus, letting ${\fancyscript{Y}}_A = \frac{\hat{p}_{D}- \zeta _0}{\zeta _1 - \zeta _0}$, we have $E[{\fancyscript{Y}}_A \mid {\mathbf{X}}] = g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}).$

Design B

In Design B, we also genotype a random sample of size $n$, but observe on everyone a perfect negative predictor $U$ satisfying $P(D=0 \mid U=0)=1.$ The EMR algorithm is developed among those individuals with $U=1.$ In addition to assumption (${\fancyscript{A}}$), we assume that $U$ is independent of ${\mathbf{X}}$ conditional on true disease status $D.$ We let ${\tilde{p}}_{D}= {\tilde{p}}_{D}(U) = \hat{p}_{D}U.$ Defining $\mu _d=E[\hat{p}_{D}\mid U=1, D=d]$ for $d=0,1$, $\rho = P(D=1 \mid U=1)$, $\pi _U = P(U=1)$, and ${\tilde{\mu }}_0 = \mu _0 \frac{\pi _U - \rho \pi _U}{1-\rho \pi _U}$, we may calculate $E[{\tilde{p}}_{D}\mid {\mathbf{X}}] = \sum _{d \in \{0,1\}} \mu _d P(U=1, D=d \mid {\mathbf{X}}) = \sum _{d \in \{0,1\}} \mu _d P(U=1 \mid D=d)P(D=d\mid {\mathbf{X}})= {\tilde{\mu }}_0\ + (\mu _1-{\tilde{\mu }}_0)g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})$ since $P(U=1 \mid D=1)=1$ and $P(U=1\mid D=0) = \frac{\pi _U - \rho \pi _U}{1-\rho \pi _U}$ by an application of Bayes rule. Thus, letting ${\fancyscript{Y}}_B = \frac{{\tilde{p}}_D - {\tilde{\mu }}_0}{\mu _1-{\tilde{\mu }}_0}$, we have $E[{\fancyscript{Y}}_B \mid {\mathbf{X}}]= g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}).$

Design C

In Design C, we first partition the full EMR into a disease-mart ($M=1$) that includes all disease cases and a control-mart ($M=0$) of disease-free individuals. We develop and apply our algorithm to calculate $\hat{p}_{D}$ only among individuals with $M=1.$ Let $\text {PPV}(S) = P(D=1 \mid M=1, \hat{p}_{D}> p_{S})$, and we assume a design with $m$ controls per case sampled from the control-mart. Let $V$ be the indicator that an individual is sampled in our study, and let ${\tilde{p}}_{D}={\tilde{p}}_{D}(M) = \hat{p}_{D}M.$

We may calculate $E[{\tilde{p}}_{D}\mid {\mathbf{X}}, D, V=1] = E[{\tilde{p}}_{D}\mid D, V=1] = DE[\hat{p}_{D}\mid D=1, M=1, \hat{p}_{D}>p_{S}]P(M=1 \mid D=1, V=1) + (1-D)E[\hat{p}_{D}\mid D=0, M=1, \hat{p}_{D}>p_{S}]P(M=1\mid D=0,V=1)= D\xi _1 + (1-D)\xi _0(1-\pi )$ where $\xi _d=E[\hat{p}_{D}\mid D=d, M=1, \hat{p}_{D}>p_{S}]$ and $\pi =P(M=0\mid D=0,V=1).$ In this calculation, we have used that ${\tilde{p}}_{D}=0$ when $M=0$ and that $P(M=1\mid D=1,V=1)=1$ because the initial partition has perfect sensitivity. We further calculate: $\pi = \frac{P(M=0, D=0 \mid V=1)}{P(D=0 \mid V=1)} = \frac{P(M=0 \mid V=1)}{P(M=0 \mid V=1) + P(D=0 \mid M=1, V=1)P(M=1 \mid V=1)} = \frac{\frac{m}{m+1}}{\frac{m}{m+1} + (1-\text {PPV}(S))\frac{1}{m+1}} = \frac{m}{m+(1-\text {PPV}(S))}.$ Then, letting ${\fancyscript{Y}}_C = \frac{{\tilde{p}}_{D}- \xi _0(1-\pi )}{\xi _1 - \xi _0(1-\pi )}$, we have that $E[{\fancyscript{Y}}_C \mid {\mathbf{X}}, V=1] = E[D \mid {\mathbf{X}}, V=1].$

Finally, using Bayes rule, we see that $E[D\mid {\mathbf{X}}, V=1] =\frac{P(V=1 \mid D=1)g(\beta ^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})}{P(V=1 \mid D=1) g(\beta ^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}) + P(V=1 \mid D=0)(1-g(\beta ^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}))} =\frac{\exp (\beta ^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})}{\lambda +\exp (\beta ^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})}$ where $\lambda = \frac{P(V=1 \mid D=0)}{P(V=1 \mid D=1)}.$ Letting $\beta _0^* = \beta _0 - \log \{\lambda \}$, we have $E[D \mid {\mathbf{X}}, V=1] = \frac{\exp (\beta _0^* + \beta _1Z +{\beta }_2^{\mathsf{\scriptscriptstyle {T}}}{\mathbf {W}})}{1+ \exp (\beta _0^* + \beta _1Z +{\beta }_2^{\mathsf{\scriptscriptstyle {T}}}{\mathbf {W}})}.$ Thus, $E[{\fancyscript{Y}}_C \mid {\mathbf{X}}, V=1] = g({{\beta }^*}^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})$, where ${\beta }^* = (\beta _0^*, \beta _1, {\beta }_2^{\mathsf{\scriptscriptstyle {T}}})^{\mathsf{\scriptscriptstyle {T}}}$.

Power and bias calculations

For simplicity we derive expressions under Design A. When using $\hat{p}_{D}$, the estimator ${\hat{{\beta }}}$ solves $U({\beta })=\frac{1}{n}\sum _{i=1}^n\psi _i({\beta })=\frac{1}{n}\sum _{i=1}^n {\mathbf{X}}_i({\fancyscript{Y}}_i - g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}_i))=0$, so $\sqrt{n}({\hat{{\beta }}}- {\beta }_0) \mathop {\rightarrow }\limits ^{{\fancyscript{D}}} N(0, V({\beta }))$ where $V({\beta })=B({\beta })^{-1}A({\beta })(B({\beta })^{-1})^{\mathsf{\scriptscriptstyle {T}}}$ where $B({\beta }) = E\left[ g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})(g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})-1) {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\right]$ and $A({\beta }) = E\left[ ({\fancyscript{Y}}- g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}))^2 {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\right] = E\left[ \left\{ \text{ Var }({\fancyscript{Y}}| D) + E[{\fancyscript{Y}}\mid D]^2 \right\} E\left[ {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\mid D \right] \right] - 2 E\left[ E\left[ {\fancyscript{Y}}\mid D \right] E\left[ g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}) {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\mid D \right] \right] + E\left[ g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})^2 {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\right]$, using assumption (${\fancyscript{A}}$). We can further expand this since ${\mathbf{X}}=(1, Z)^{\mathsf{\scriptscriptstyle {T}}}$, for SNP $Z;$ in particular, for any function $f$, $E[f(Z) \mid D] = \frac{D}{P(D=1)} E[f(Z)g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})] + \frac{1-D}{P(D=0)} E[f(Z)(1-g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}))].$ Letting $\mu _d=E[{\fancyscript{Y}}\mid D=d]$ and $\xi _d=\text{ Var }({\fancyscript{Y}}\mid D=d)$, we can rewrite $A({\beta })$ as: $A({\beta }) = \left\{ \xi _0 + \mu _0^2 \right\} E[{\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}] + \left\{ \xi _1 + \mu _1^2 - \xi _0 -\mu _0^2 -2 \mu _0 \right\} E[g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}){\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}] + (2(\mu _0 - \mu _1) +1) E[g({\beta }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})^2{\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}]$.

To compare power to results using ${\tilde{D}}$ in the misspecified model, we now consider the distribution of ${\hat{{\gamma }}}$ which solves $U({\gamma })=\frac{1}{n}\sum _{i=1}^n\psi _i({\gamma })=\frac{1}{n}\sum _{i=1}^n {\mathbf{X}}_i({\tilde{D}}_i - g({\gamma }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}_i))=0.$ Then $\sqrt{n}({\hat{{\gamma }}}- {\gamma }^*({\beta })) \rightarrow N(0, V^*({\beta }))$, where ${\gamma }^*({\beta })$ is a constant. To estimate ${\gamma }^*$, we can proceed as in Neuhaus (1999) and use results from work on misspecified models to see that estimates from the false model $P_F({\tilde{D}}=1 \mid Z) = g(\gamma _0 + \gamma _1Z)$ converge to values $(\gamma _0^*, \gamma _1^*)$ that minimize the Kullback–Leibler divergence between the false model and the true model $P_T({\tilde{D}}=1 \mid Z) = (1-S) + (\text{ SE }(S) + S - 1)g(\beta _0 + \beta _1Z)$ (Neuhaus 1999; Kullback 1959). The Kullback–Leibler divergence between these two models is $E_Z[ \log \{ \frac{P_T({\tilde{D}}=1 \mid Z)}{P_F({\tilde{D}}=1 \mid Z)} \}P_T({\tilde{D}}=1 \mid Z)+ \log \{ \frac{P_T({\tilde{D}}=0 \mid Z)}{P_F({\tilde{D}}=0 \mid Z)} \} P_T({\tilde{D}}=0 \mid Z)].$ By taking derivatives with respect to $\gamma _0$ and $\gamma _1$ and setting them to 0, we find two equations: $0= \left\{ \alpha _0 + \alpha _1g(\beta _0) -g(\gamma _0)\right\} (1-p_Z)^2 + \left\{ \alpha _0 + \alpha _1g(\beta _0+\beta _1) -g(\gamma _0+\gamma _1)\right\} 2p_Z(1-p_Z) + \left\{ \alpha _0 + \alpha _1g(\beta _0+2\beta _1) -g(\gamma _0+2\gamma _1)\right\} p_Z^2$ and $0= \left\{ \alpha _0 + \alpha _1g(\beta _0+\beta _1)\right.$ $\left. -g(\gamma _0+\gamma _1)\right\} 2p_Z(1-p_Z) +2 \left\{ \alpha _0 + \alpha _1g(\beta _0+2\beta _1) -g(\gamma _0+2\gamma _1)\right\} p_Z^2$, where $\alpha _0=1-S$, $\alpha _1=\text{ SE }(S)+S-1.$ Here, we assume that the SNP $Z\sim \text {Bin}(2, p_Z)$, where $p_Z$ is the MAF. Simultaneously solving these two equations for $(\gamma _0, \gamma _1)$ yields the desired $(\gamma ^*_0, \gamma ^*_1).$ The calculation of $V^*({\gamma }^*)$ proceeds similarly to the calculation of $V({\beta }).$ $V^*({\gamma }^*)=B^*({\gamma }^*)^{-1}A^*({\gamma }^*)(B^*({\gamma }^*)^{-1})^{\mathsf{\scriptscriptstyle {T}}}.$ where $B^*({\gamma }) = E\left[ g({\gamma }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})(g({\gamma }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})-1) {\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}\right]$ as before. Here, though, $A^*({\gamma })= (1-S)E[{\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}]+(\text{ SE }(S)-3(1-S)) E[g({\gamma }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}}){\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}] +(2(1-S-\text{ SE }(S)) + 1) E[g({\gamma }^{\mathsf{\scriptscriptstyle {T}}}{\mathbf{X}})^2{\mathbf{X}}{\mathbf{X}}^{\mathsf{\scriptscriptstyle {T}}}].$