An Empirical Bayes risk prediction model using multiple traits for sequencing data

Appendix

A. the best linear unbiased predictor

1. q is a standardized vector, and E(q)=0 and Var(q)=1.

2. ε_i have mean zero and variance ψ_i. Normally, ψ_i=1.

3. assuming a simple linear regression model (Equation 5), we have,

   E(Ci)=E(βi,vq+εi)=0Var(Ci)=Var(βi,vq+εi)=βi,v2+ψi

   Extending C_i to C, we get

   E(C)=0Var(C)=βvβ′v+Ψ

   where Ψ is the variance-covariance matrix of ε.

4. The joint distribution of (q, C) is (qC)~[(00), (1βvβvβvβ′v+Ψ)] since,

   Cov(q, C)=Cov(q, βvq+ε)=Cov(q, βvq)+Cov(q, ε)=βvCov(q, q)+0=βvCor(q, C)=βv1×(βvβ′v+Ψ)1/2

   where q and ε are uncorrelated.

5. The best linear predictor is given by

   g(X)=μy+ρσyσx(X−μx)

The predicted value of q is denoted as q^# based on g(X). q and C correspond to y and x in the above equation. According to (i) and (iii),

E(q)=0⇒μy=0E(C)=0⇒μx=0Var(q)=1⇒σy=1Var(C)=βvβv′+Ψ⇒σx=(βvβv′+Ψ)1/2Cor(q, C)=βv1×(βvβv′+Ψ)1/2⇒ρ=βv(βvβv′+Ψ)1/2

Finally, the best predictor equation becomes,

q#=ρσyσxX=β′v(βvβ′v+Ψ)−1C=11+λv2β′vΨ−1C

where Ψ is the variance-covariance matrix of ε, and λv2=β′vΨ−1βv which is a squared multivariate distance of β_v from the origin.

B: thetstatistic

If ε_i follows a normal distribution, the t-test statistic for the null hypothesis β_i,v=0 under model (5) follows a non-central t distribution,

Ti,v=β^i,v−0se(β^i,v)=2c0β^i,v∼tn−2(γi,v) where c02=∑j=1n(qj−q¯)24

Proof: Given q, the ordinary least squares estimate of coefficient β_i,v is β^i,v|q=(q′q)−1q′Ci (Montgomery et al., 2006),

Var(β^i,v|q)=(q′q)−1q′q(q′q)−1Var(Ci|q)=(q′q)−1se(β^i,v|q)=(q′q)−1/2

It is known that E(q)=0,

q′q=∑j=1n(qj−q¯)2

If ε_i is normally distributed, the t-test statistic follows,

Ti,v=β^i,v−0se(β^i,v)=(q′q)1/2β^i,v=2c0β^i,v∼tn−2(γi,t)

where c02=∑j=1n(qj−q¯)24, then 2c0=(∑j=1n(qj−q¯)2)1/2=(q′q)1/2.

C: proof of proposition 1

The proof follows the results in literature (Ker, 2001; Nadarajah and Kotz, 2008). To simplify the notation, we use Z₁ and Z₂ to represent Z_i,d and Z_i,v, respectively. Assume that (Z₁, Z₂) follow a bivariate normal distribution (Z1Z2)∼MVN(γ, Σ) where γ=(γ1, γ2)T, Σ=(σ12pσ1σ2pσ1σ2σ22) where ρ denotes the correlation between Z₁ and Z₂. Let Z=min(Z₁, Z₂) and I=r if Z=Z_r (r=1, 2).

Lemma 1: Let ϕ(.) and Φ(.) are the univariate standard normal density (pdf) and its corresponding cumulative distribution function (cdf), respectively. Then the conditional pdf of Z given I=r for all z and r takes 1 or 2 is,

fr(z)=Pr(Z≤z|I=r){pr−1ϕ(z−γrσr)(1−Φ(z−γr∗σr∗))if  ρ≠1.ϕ(z−γrσr)otherwise.

where γr∗=αrγ3−r+(1−αr)γr, σr∗=αrσ3−r1−ρ2, and αr=11−ρσ3−rσr. Here, Z=min(Z₁, Z₂) works for ρ≠1. If ρ=1, Z=Z₁=Z₂.

And pr=P(I=r)=Φ(γr−γ3−rσ12+σ22−2ρσ1σ2)

Proof: Let ℵ be the distribution of pair (Z, I), then

ℵ(z, 1)=Pr(Z≤z, I=1)=Pr(I=1)−Pr(Z2>Z1>z)=p1−∫z∞∫z1∞ϕ((z2−γ2)−ρσ2σ1(z1−γ1)σ21−ρ2)ϕ(z1−γ1σ1)dz2dz1=p1−∫z∞(1−Φ((z−γ2)−ρσ2σ1(z−γ1)σ21−ρ2))ϕ(z−γ1σ1)dz

Hence, for all z,

ℵ(z, 1)=Pr(Z≤z, I=1)=Pr(Z≤z|I=1)Pr(I=1)=p1f1(z)=1σ1ϕ(z−γ1σ1)(1−Φ((z−γ2)−ρσ2σ1(z−γ1σ1)σ21−ρ2))=1σ1ϕ(z−γ1σ1)Φ(ρ(z−γ1)σ11−ρ2−z−γ2σ21−ρ2)

where Φ(a)=1–Φ(a).

The equation of p₁ is derived from,

p1=Pr(Z1−Z2<0)=∫−∞γ1−γ2σ12+σ22−2ρσ1σ2ϕ(z)dz

Proposition 1: The distribution of the pair (Z, I) uniquely determines the distribution of the pair (Z₁, Z₂).

Proof: Suppose (Z, I) has the distribution ℵ, and we use parameters (γ′1, γ′2) and σ′r for r=1 or 2 to denote any bivariate normal pdf of (Z, I). According to the above results,

limz→−∞prfr(z)ϕ(z−γrσr)=limz→−∞prfr(z)ϕ(z−γ′rσ′r)=1

Hence, ϕ(z−γrσr)ϕ(z−γ′rσ′r)⇒1, then we will get γr=γ′r, σr=σ′r and ρ=ρ′.

Similarly, for I=2, the pdf of (Z, I) is derived by

ℵ(z, 2)=Pr(Z≤z, I=2)=p2f2(z)=1σ2ϕ(z−γ2σ2)Φ(ρ(z−γ2)σ21−ρ2−z−γ1σ11−ρ2)

Finally, the probability density function of Z=min(Z₁, Z₂) for I=1 or 2 will be given by,

ℵ(z, 1)+ℵ(z, 2)=p1f1(z)+p2f2(z)=1σ1ϕ(z−γ1σ1)Φ(ρ(z−γ1)σ11−ρ2−z−γ2σ21−ρ2)+1σ2ϕ(z−γ2σ2)Φ(ρ(z−γ2)σ21−ρ2−z−γ1σ11−ρ2)

when Z₁ and Z₂ are replaced by Z_i,d and Z_i,v separately, γ₁ and γ₂ denote γ_i,d and γ_i,v, and Z is replaced by Z_i,m, The pdf of Z_i,m=min(Z_i,d, Z_i,v) is derived.

In the same way, we can derive the pdf of Z=max(Z₃, Z₄). Assume that (Z₃, Z₄) follow a bivariate normal distribution, (Z3Z4)∼MVN(γ, Σ) where γi=(γ3, γ4)T, Σ=(σ32ρσ3σ4ρσ3σ4σ42) where ρ denotes the correlation between Z₃ and Z₄. Let Z=max(Z₃, Z₄) and I=r if Z=Z_r (r=3 or 4).

Lemma 2: Let ϕ(.) and Φ(.) are the univariate standard normal density (pdf) and its corresponding cumulative distribution function (cdf), respectively. Then the pdf of Z given I=r for all z and r takes 3 or 4 is,

fr(z)=Pr(Z≤z|I=r){pr−1ϕ(γr−zσr)Φ(z−γr∗σr∗)if  ρ≠1.ϕ(γr−zσr),otherwise.

where γr∗=αrγ7−r+(1−αr)γr, σr∗=αrσ7−r1−ρ2, and αr=11−ρσ7−rσr.

Proof: let ℵ be the distribution of pair (Z, I), then

ℵ(z, 3)=Pr(Z≤z, I=3)=Pr(Z4<Z3≤z)=∫−∞z∫−∞z3ϕ((z3−γ3)−ρσ3σ4(z3−γ4)σ31−ρ2)ϕ(z4−γ4σ4)dz4dz3=∫−∞zΦ(ρσ3σ4(γ4−z)−(γ3−z)σ31−ρ2)ϕ(z−γ4σ4)dz

Then the pdf of (Z, 3) is given by

p3f3(z)=1σ4ϕ(γ4−zσ4)Φ(ρ(γ4−z)σ41−ρ2−γ3−zσ31−ρ2)

where ϕ(a)=ϕ(–a).

Similarly, the pdf of (Z, 4) is given by,

p4f4(z)=1σ3ϕ(γ3−zσ3)Φ(ρ(γ3−z)σ31−ρ2−γ4−zσ41−ρ2)

The equation of p₃ is derived from,

p3=Pr(Z3−Z4>0)=∫γ4−γ3σ32+σ42−2ρσ3σ4)∞ϕ(z)dz

Proposition 2: The distribution of the pair (Z, I) uniquely determines the distribution of the pair (Z₃, Z₄).

Proof: same as the above proposition 1.

Finally, the pdf of Z=max(Z₁, Z₂) is given by,

p3f3(z)+p4f4(z)=1σ4ϕ(γ4−zσ4)Φ(ρ(γ4−z)σ41−ρ2−γ3−zσ31−ρ2)+1σ3ϕ(γ3−zσ3)Φ(ρ(γ3−z)σ31−ρ2−γ4−zσ41−ρ2)

Again, when Z₃ and Z₄ are replaced by Z_i,v and Z_i,d separately, γ₃ and γ₄ denote γ_i,v and γ_i,d, and Z is replaced by Z_i,m, The pdf of Z_i,m=max(Z_i,d, Z_i,v) is derived.

D: Effects of multiple quantitative traits

GAW 17 provides three quantitative traits, such as: Q1, Q2 and Q4. In particular, the residual heritability of Q1 is 0.44, and covariates, age and smoking status, act on the Q1. Q2 has a heritability of 0.29, but is not affected by covariates. Q4 is simulated by various protective genes, and it has a heritability of 0.7. According to the disease risk simulation, Q1 and Q2 are positively related to the disease trait, but Q4 is negatively correlated to the disease status (Appendix Table A1). When we add three quantitative traits in the disease risk EB classifier separately (Appendix Table A2), five true genes are detected by adding Q1 in the EB model, four true genes identified through incorporating Q2 in the EB model, and only two true genes explored by adding Q4 in the EB method. In particular, all three quantitative traits detect FLT1 gene, but Q1 provides the highest ranking position compared to Q2 and Q4. Appendix Table A3 further evaluates the performances of three quantitative traits on predicting disease risk. It clearly illustrates that adding Q1 in the EB model gives the smallest cross-validation error (0.26 with SD 0.0003) and the largest AUC (0.91 with SD 0.022), followed by Q2, then Q4. The results agree with the simulation setting for each quantitative trait.

E: Multiple replicates

The GAW 17 sequence data contains R replicated phenotypes. This results in the analogous Z statistic for every replicate. According to the Brown and Stein Bayesian model, each Z_i,a (a=1, …, R) is derived from a normal distribution N(γi, σi2), then the expectation of a sequence of statistics (Z̅_i) is still γ_i, that is E(Z¯i)=E(Zi,1+…+Zi,RR)=γi. But the variance of this new estimator deviates from the unit variance due to the correlations among the sequence of Z statistics, which may violate the assumption of the statistic before estimating EB estimates. To interpret the variance of Z̅_i, we define a correlation r for any two replicates of statistics (Z_i,a and Z_i,b) where a, b=1, …, R and a≠b, then Cor(Z_i,a, Z_i,b)=r.

Var(Z¯i)=Var(Zi,1+…+Zi,RR)=Var(Zi,1)+…+Var(Zi,R)+(R×(R−1))rR2=1+(R−1)rRσi2

It is clearly seen that the variance of Z̅_i is not a constant, but the parameter correlation r can be empirically estimated by

Zi,a∼N(γi, σi2)   and   Z¯i∼N(γi, 1+(R−1)rRσi2)Cor(Zi,a, Zi,b)=r   for any a, b, i

where a and b measure different replicates, and i indicates the ith feature. Because σi2 is close to 1, we could let σi2 take 1 to simplify our calculation without lose the accuracy of estimator.

A new statistic, Zi‡=Z¯i1+(R−1)rR, defined from Z̅_i has a unite variance, and follows a normal distribution with mean γi‡=γi1+(R−1)rR. The relevant Bayesian estimate of parameter γi‡ is achieved by E(γi‡|Zi‡) via shrinking the estimator Zi‡. The Empirical Bayes estimate of true γ_i is estimated by E(γi|Zi‡)=E(γi‡|Zi‡)1+(R−1)rR where E(γi‡|Zi‡)=E(γi1+(R−1)rR|Zi‡)=E(γi|Zi‡)1+(R−1)rR. Then we incorporate estimates E(γi|Zi‡) into the prediction model. In particular, r is estimated by

∑i=1N∑a=1R(Zi,a−Z¯i)2R−1N=1−rr^=1−∑i=1NsZi2N

where sZi2 is the sample variance of a sequence of statistics {Z_i,1, …, Z_i,R}at the ith feature, and i=1, …, N.