From 026fa1a09ff5aca4776900bcb9d7f45c2d94ca86 Mon Sep 17 00:00:00 2001 From: nicholst Date: Wed, 23 Mar 2016 10:47:47 +0000 Subject: [PATCH] All edits delivered by Xu on 2015/12/30. Code now accepts AE/ACE option. Various other clean ups. --- apace/ACEfit.m | 286 ++++++++++++++++++++----------------- apace/ACEfit_Perm.m | 198 +++++++++++++------------ apace/ACEfit_Results.m | 198 ++++++++++--------------- apace/APACEsummary.m | 8 +- apace/AgHe_Method.m | 101 ++++++------- apace/Boot_CIs.m | 143 ++++++++++--------- apace/CheckData.m | 3 +- apace/PrepParallel.m | 20 +-- apace/ProcData.m | 14 +- apace/Reorder.m | 53 ++++++- apace/WK_ParallelOutline.m | 218 +++++++++++++++------------- 11 files changed, 651 insertions(+), 591 deletions(-) diff --git a/apace/ACEfit.m b/apace/ACEfit.m index b5884fa..e1d0720 100644 --- a/apace/ACEfit.m +++ b/apace/ACEfit.m @@ -5,9 +5,9 @@ % % The LRT statistic is chosen as the test statistic rather than h2 % estimator as theoretically h2 estimator is not sufficient nor pivotal -% while LRT statistic is pivotal. The distribution of LRT statistic does not -% depend on unknown nuisance parameters and remains the same under the null -% hypothesis, and thus the pivotal test statistic is recommendated in +% while LRT statistic is pivotal. The distribution of LRT statistic does +% not depend on unknown nuisance parameters and remains the same under the +% null hypothesis, and thus the pivotal test statistic is recommendated in % permutation test. % @@ -29,7 +29,7 @@ Sigma2 = sum(Res.^2,2)/(n-size(X,2)); [XTX,pinvXTX] = Mk_XTXpinv(size(I_MZ,1),size(I_DZ,1),size(I_Sib,1),n); - +BadElm = 0; switch upper(ACEfit_Par.Model) @@ -47,7 +47,7 @@ iY = Y(i,:)'; iRes = Res(i,:)'; sig2 = Sigma2(i); - + % % Heritability estimation % @@ -63,91 +63,102 @@ % % Likelihood ratio test % - if ( min(HH(:,3))>=0 ) - hh0 = HH(:,3); - else - hh0 = HH(:,1); - end - if ( min(HH(:,2))>=0 ) - hh_AE = HH(:,2); - else - hh_AE = HH(:,1); - end - hh0 = hh0/sum(hh0); - hh_AE = hh_AE/sum(hh_AE); - [det_Va,det_V0,det_V_AE] = deal(1); - [invMZR,invOTR,invR0,invMZR_AE,invOTR_AE] = deal(cell(max(nFam),1)); - [invMZR{1},invOTR{1},invR0{1},invMZR_AE{1},invOTR_AE{1}] = deal(1); - - for nF=2:max(nFam) - n_nF = sum(nFam==nF); - n_nFMZ = sum(nFam(1:nMZF)==nF); + if (hha(3)<1e-5) + % Zero environmental variance... give up + BadElm = BadElm + 1; + AsyPval_A(i) = 1; + AsyPval_C(i) = 1; + Stats(i) = 0; + Stats_C(i) = 0; + else + + if ( min(HH(:,3))>=0 ) + hh0 = HH(:,3); + else + hh0 = HH(:,1); + end + if ( min(HH(:,2))>=0 ) + hh_AE = HH(:,2); + else + hh_AE = HH(:,1); + end + hh0 = hh0/sum(hh0); + hh_AE = hh_AE/sum(hh_AE); + + [det_Va,det_V0,det_V_AE] = deal(1); + [invMZR,invOTR,invR0,invMZR_AE,invOTR_AE] = deal(cell(max(nFam),1)); + [invMZR{1},invOTR{1},invR0{1},invMZR_AE{1},invOTR_AE{1}] = deal(1); + + for nF=2:max(nFam) + n_nF = sum(nFam==nF); + n_nFMZ = sum(nFam(1:nMZF)==nF); + + %%% ACE model + a = hha(1)/2+hha(3); + b = hha(1)/2+hha(2); + R = eye(nF)*a + ones(nF)*b; + det_Va = det_Va*det(R)^(n_nF-n_nFMZ); + invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + + R(1,2) = hha(1)+hha(2); + R(2,1) = R(1,2); + det_Va = det_Va*det(R)^n_nFMZ; + invMZR{nF} = sparse(pinv(R)); + + %%% CE model + R0 = eye(nF)*hh0(3) + ones(nF)*hh0(2); + det_V0 = det_V0*det(R0)^n_nF; + invR0{nF} = (eye(nF) - ones(nF)*hh0(2)/(hh0(3)+hh0(2)*nF))/hh0(3); + + %%% AE model + a = hh_AE(1)/2+hh_AE(3); + b = hh_AE(1)/2; + R_AE = eye(nF)*a + ones(nF)*b; + det_V_AE = det_V_AE*det(R_AE)^(n_nF-n_nFMZ); + invOTR_AE{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + + R_AE(1,2) = hh_AE(1); + R_AE(2,1) = R_AE(1,2); + det_V_AE = det_V_AE*det(R_AE)^n_nFMZ; + invMZR_AE{nF} = sparse(pinv(R_AE)); + end - %%% ACE model - a = hha(1)/2+hha(3); - b = hha(1)/2+hha(2); - R = eye(nF)*a + ones(nF)*b; - det_Va = det_Va*det(R)^(n_nF-n_nFMZ); - invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + % ACE model + Ua = blkdiag(invMZR{nFam(1:nMZF)},invOTR{nFam(nMZF+1:end)}); + Za = pinv(X'*Ua*X); + Pa = Ua-Ua*X*Za*X'*Ua; + rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; - R(1,2) = hha(1)+hha(2); - R(2,1) = R(1,2); - det_Va = det_Va*det(R)^n_nFMZ; - invMZR{nF} = sparse(pinv(R)); + % CE model + U0 = blkdiag(invR0{nFam}); + Z0 = pinv(X'*U0*X); + P0 = U0-U0*X*Z0*X'*U0; + rl0 = -(iY'*P0*iY/sig2+log(det_V0)+log(det(X'*U0*X)))/2; - %%% CE model - R0 = eye(nF)*hh0(3) + ones(nF)*hh0(2); - det_V0 = det_V0*det(R0)^n_nF; - invR0{nF} = (eye(nF) - ones(nF)*hh0(2)/(hh0(3)+hh0(2)*nF))/hh0(3); + % AE model + U_AE = blkdiag(invMZR_AE{nFam(1:nMZF)}, invOTR_AE{nFam(nMZF+1:end)}); + Z_AE = pinv(X'*U_AE*X); + P_AE = U_AE-U_AE*X*Z_AE*X'*U_AE; + rl_AE = -(iY'*P_AE*iY/sig2+log(det_V_AE)+log(det(X'*U_AE*X)))/2; - %%% AE model - a = hh_AE(1)/2+hh_AE(3); - b = hh_AE(1)/2; - R_AE = eye(nF)*a + ones(nF)*b; - det_V_AE = det_V_AE*det(R_AE)^(n_nF-n_nFMZ); - invOTR_AE{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + % Compute LRT test statistic + T = -2*(rl0-rla); + T_C = -2*(rl_AE-rla); + + % Obtain the resulting asymptotic p-value + if T>0 + AsyPval_A(i) = (1-chi2cdf(T,1))/2; + end + if T_C>0 + AsyPval_C(i) = (1-chi2cdf(T_C,1))/2; + end + + % Save LRT test statistic + Stats(i) = max(0,T); + Stats_C(i) = max(0,T_C); - R_AE(1,2) = hh_AE(1); - R_AE(2,1) = R_AE(1,2); - det_V_AE = det_V_AE*det(R_AE)^n_nFMZ; - invMZR_AE{nF} = sparse(pinv(R_AE)); - end - - % ACE model - Ua = blkdiag(invMZR{nFam(1:nMZF)},invOTR{nFam(nMZF+1:end)}); - Za = pinv(X'*Ua*X); - Pa = Ua-Ua*X*Za*X'*Ua; - rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; - - % CE model - U0 = blkdiag(invR0{nFam}); - Z0 = pinv(X'*U0*X); - P0 = U0-U0*X*Z0*X'*U0; - rl0 = -(iY'*P0*iY/sig2+log(det_V0)+log(det(X'*U0*X)))/2; - - % AE model - U_AE = blkdiag(invMZR_AE{nFam(1:nMZF)}, invOTR_AE{nFam(nMZF+1:end)}); - Z_AE = pinv(X'*U_AE*X); - P_AE = U_AE-U_AE*X*Z_AE*X'*U_AE; - rl_AE = -(iY'*P_AE*iY/sig2+log(det_V_AE)+log(det(X'*U_AE*X)))/2; - - % Compute LRT test statistic - T = -2*(rl0-rla); - T_C = -2*(rl_AE-rla); - - % Obtain the resulting asymptotic p-value - if T>0 - AsyPval_A(i) = (1-chi2cdf(T,1))/2; - end - if T_C>0 - AsyPval_C(i) = (1-chi2cdf(T_C,1))/2; end - - % Save LRT test statistic - Stats(i) = max(0,T); - Stats_C(i) = max(0,T_C); - end end @@ -211,11 +222,11 @@ % if isnan(mGq3e2O) % mGq3e2O = mean(e2O(e2O>=q3e2O)); % end - + % Summary statistics: meane2, we2, median, q3, mean(e2>median), mean(e2>q3) ACEfit_Par.SummaryE = [meane2O; we2O; mede2O; q3e2O; mGmede2O; mGq3e2O]; - + % Write out the ouptut images WriteData(h2, ACEfit_Par.Vs, 'ACE_A_h2', ACEfit_Par.ResDir); WriteData(c2, ACEfit_Par.Vs, 'ACE_C_c2', ACEfit_Par.ResDir); @@ -239,7 +250,7 @@ iY = Y(i,:)'; iRes = Res(i,:)'; sig2 = Sigma2(i); - + % % Heritability estimation % @@ -249,63 +260,72 @@ else hha = HH(:,1); end - hha = hha/sum(hha); + hha = hha/sum(hha); - h2(i) = hha(1); - e2(i) = hha(3); + h2(i) = hha(1); + e2(i) = hha(3); if ~ACEfit_Par.NoImg % % Likelihood ratio test % - index = ACEcode(hha); - if ( index==ACEcode([0 0 1]) ) - % Parameter A is zero, test-statistic is zero - T = 0; - else - det_Va = 1; - [invMZR,invOTR] = deal(cell(max(nFam),1)); - [invMZR{1},invOTR{1}] = deal(1); + + if (hha(3)<1e-5) + % Zero enviornmental variance... give up + BadElm = BadElm + 1; + AsyPval_A(i) = 1; + Stats(i) = 0; + else - for nF=2:max(nFam) - n_nF = sum(nFam==nF); - n_nFMZ = sum(nFam(1:nMZF)==nF); + index = ACEcode(hha); + if ( index==ACEcode([0 0 1]) ) + % Parameter A is zero, test-statistic is zero + T = 0; + else + det_Va = 1; + [invMZR,invOTR] = deal(cell(max(nFam),1)); + [invMZR{1},invOTR{1}] = deal(1); - %%% AE model - a = hha(1)/2+hha(3); - b = hha(1)/2; - Ra = eye(nF)*a + ones(nF)*b; - det_Va = det_Va*det(Ra)^(n_nF-n_nFMZ); - invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + for nF=2:max(nFam) + n_nF = sum(nFam==nF); + n_nFMZ = sum(nFam(1:nMZF)==nF); + + %%% AE model + a = hha(1)/2+hha(3); + b = hha(1)/2; + Ra = eye(nF)*a + ones(nF)*b; + det_Va = det_Va*det(Ra)^(n_nF-n_nFMZ); + invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + + Ra(1,2) = hha(1); + Ra(2,1) = Ra(1,2); + det_Va = det_Va*det(Ra)^n_nFMZ; + invMZR{nF} = sparse(pinv(Ra)); + end - Ra(1,2) = hha(1); - Ra(2,1) = Ra(1,2); - det_Va = det_Va*det(Ra)^n_nFMZ; - invMZR{nF} = sparse(pinv(Ra)); + % AE model + Ua = blkdiag(invMZR{nFam(1:nMZF)}, invOTR{nFam(nMZF+1:end)}); + Za = pinv(X'*Ua*X); + Pa = Ua-Ua*X*Za*X'*Ua; + rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; + + % E model + rl0 = -(iY'*P0*iY/sig2+log(det(X'*X)))/2; + + % Compute LRT test statistic + T = -2*(rl0-rla); end - % AE model - Ua = blkdiag(invMZR{nFam(1:nMZF)}, invOTR{nFam(nMZF+1:end)}); - Za = pinv(X'*Ua*X); - Pa = Ua-Ua*X*Za*X'*Ua; - rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; + % Obtain the resulting asymptotic p-value + if T>0 + AsyPval_A(i) = (1-chi2cdf(T,1))/2; + end + + % Save LRT test statistic + Stats(i) = max(0,T); - % E model - rl0 = -(iY'*P0*iY/sig2+log(det(X'*X)))/2; - - % Compute LRT test statistic - T = -2*(rl0-rla); - end - - % Obtain the resulting asymptotic p-value - if T>0 - AsyPval_A(i) = (1-chi2cdf(T,1))/2; end - - % Save LRT test statistic - Stats(i) = max(0,T); - end end @@ -359,6 +379,10 @@ WriteData(sqrt(Sigma2), ACEfit_Par.Vs, 'Stdev', ACEfit_Par.ResDir); end +if BadElm>0 + fprintf('WARNING: Voxels/elements found with zero environmental variance (%d). Check for serverely sparse discrete data.\n',BadElm); +end + % Plot the distribution of overall h2 estimates f = 0; % 1) h2_dist @@ -398,13 +422,13 @@ WriteData(Stats, ACEfit_Par.Vs, 'ACE_A_LRT', ACEfit_Par.ResDir); WriteData(Stats_C, ACEfit_Par.Vs, 'ACE_C_LRT', ACEfit_Par.ResDir); WriteData(-log10(AsyPval_A), ACEfit_Par.Vs, 'ACE_A_LRT_vox_Pasym', ACEfit_Par.ResDir); - WriteData(-log10(AsyPval_C), ACEfit_Par.Vs, 'ACE_C_LRT_vox_Pasym', ACEfit_Par.ResDir); + WriteData(-log10(AsyPval_C), ACEfit_Par.Vs, 'ACE_C_LRT_vox_Pasym', ACEfit_Par.ResDir); case 'AE' % Save maximum LRT statistic (mT) ACEfit_Par.mT = max(Stats); ACEfit_Par.Stats = Stats; WriteData(Stats, ACEfit_Par.Vs, 'ACE_A_LRT', ACEfit_Par.ResDir); - WriteData(-log10(AsyPval_A), ACEfit_Par.Vs, 'ACE_A_LRT_vox_Pasym', ACEfit_Par.ResDir); + WriteData(-log10(AsyPval_A), ACEfit_Par.Vs, 'ACE_A_LRT_vox_Pasym', ACEfit_Par.ResDir); end diff --git a/apace/ACEfit_Perm.m b/apace/ACEfit_Perm.m index 88559c1..e134a45 100644 --- a/apace/ACEfit_Perm.m +++ b/apace/ACEfit_Perm.m @@ -37,6 +37,8 @@ Stats = zeros(nElm,1); end +BadElm = 0; + switch upper(ACEfit_Par.Model) case 'ACE' @@ -50,7 +52,7 @@ % % Heritability estimation % - [hha,HH] = LRSD_ACE_FAM(I_MZ,I_DZ,I_Sib,iRes,XTX,pinvXTX,sig2); + [hha,HH] = LRSD_ACE_FAM(I_MZ,I_DZ,I_Sib,iRes,XTX,pinvXTX,sig2); hha = hha/sum(hha); h2(i) = hha(1); @@ -60,68 +62,76 @@ % % Likelihood ratio test % - index = ACEcode(hha); - if ( index==ACEcode([0 0 1]) || index==ACEcode([0 1 1]) ) - % Parameter A is zero, test-statistic is zero - T = 0; + + if (hha(3)<1e-5) + % Zero environmental variance... give up + BadElm = BadElm + 1; + Stats(i) = 0; else - % Parameter A is non-zero... compute log-likelihood... - if ( min(HH(:,3))>=0 ) - hh0 = HH(:,3); - else - hh0 = HH(:,1); - end - hh0 = hh0/sum(hh0); - - [det_Va,det_V0] = deal(1); - [invMZR,invOTR,invR0] = deal(cell(max(nFam),1)); - [invMZR{1},invOTR{1},invR0{1}] = deal(1); - for nF=2:max(nFam) - n_nF = sum(nFam==nF); - n_nFMZ = sum(nFam(Perm_label(1:nMZF))==nF); + index = ACEcode(hha); + if ( index==ACEcode([0 0 1]) || index==ACEcode([0 1 1]) ) + % Parameter A is zero, test-statistic is zero + T = 0; + else + % Parameter A is non-zero... compute log-likelihood... + if ( min(HH(:,3))>=0 ) + hh0 = HH(:,3); + else + hh0 = HH(:,1); + end + hh0 = hh0/sum(hh0); - %%% ACE model - a = hha(1)/2+hha(3); - b = hha(1)/2+hha(2); - R = eye(nF)*a + ones(nF)*b; - det_Va = det_Va*det(R)^(n_nF-n_nFMZ); - invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + [det_Va,det_V0] = deal(1); + [invMZR,invOTR,invR0] = deal(cell(max(nFam),1)); + [invMZR{1},invOTR{1},invR0{1}] = deal(1); - R(1,2) = hha(1)+hha(2); - R(2,1) = R(1,2); - det_Va = det_Va*det(R)^n_nFMZ; - invMZR{nF} = sparse(pinv(R)); + for nF=2:max(nFam) + n_nF = sum(nFam==nF); + n_nFMZ = sum(nFam(Perm_label(1:nMZF))==nF); + + %%% ACE model + a = hha(1)/2+hha(3); + b = hha(1)/2+hha(2); + R = eye(nF)*a + ones(nF)*b; + det_Va = det_Va*det(R)^(n_nF-n_nFMZ); + invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + + R(1,2) = hha(1)+hha(2); + R(2,1) = R(1,2); + det_Va = det_Va*det(R)^n_nFMZ; + invMZR{nF} = sparse(pinv(R)); + + %%% CE model + R0 = eye(nF)*hh0(3) + ones(nF)*hh0(2); + det_V0 = det_V0*det(R0)^n_nF; + invR0{nF} = (eye(nF) - ones(nF)*hh0(2)/(hh0(3)+hh0(2)*nF))/hh0(3); + end - %%% CE model - R0 = eye(nF)*hh0(3) + ones(nF)*hh0(2); - det_V0 = det_V0*det(R0)^n_nF; - invR0{nF} = (eye(nF) - ones(nF)*hh0(2)/(hh0(3)+hh0(2)*nF))/hh0(3); - end - - % ACE model - Ua = blkdiag(invOTR{nFam}); - for nF = Perm_label(1:nMZF) - IndF = nFam(nF); - i0 = sum(nFam(1:nF-1)); - Ua(i0+[1:IndF], i0+[1:IndF]) = invMZR{IndF}; + % ACE model + Ua = blkdiag(invOTR{nFam}); + for nF = Perm_label(1:nMZF) + IndF = nFam(nF); + i0 = sum(nFam(1:nF-1)); + Ua(i0+[1:IndF], i0+[1:IndF]) = invMZR{IndF}; + end + Za = pinv(X'*Ua*X); + Pa = Ua-Ua*X*Za*X'*Ua; + rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; + + % CE model + U0 = blkdiag(invR0{nFam}); + Z0 = pinv(X'*U0*X); + P0 = U0-U0*X*Z0*X'*U0; + rl0 = -(iY'*P0*iY/sig2+log(det_V0)+log(det(X'*U0*X)))/2; + + % Compute LRT statistic + T = max(0,-2*(rl0-rla)); end - Za = pinv(X'*Ua*X); - Pa = Ua-Ua*X*Za*X'*Ua; - rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; - % CE model - U0 = blkdiag(invR0{nFam}); - Z0 = pinv(X'*U0*X); - P0 = U0-U0*X*Z0*X'*U0; - rl0 = -(iY'*P0*iY/sig2+log(det_V0)+log(det(X'*U0*X)))/2; + Stats(i) = T; - % Compute LRT statistic - T = max(0,-2*(rl0-rla)); end - - Stats(i) = T; - end end @@ -153,48 +163,56 @@ % % Likelihood ratio test % - index = ACEcode(hha); - if ( index==ACEcode([0 0 1]) ) - % Parameter A is zero, test-statistic is zero - T = 0; + + if (hha(3)<1e-5) + % Zero environmental variance... give up + BadElm = BadElm + 1; + Stats(i) = 0; else - % Parameter A is non-zero... compute log-likelihood... - det_Va = 1; - [invMZR,invOTR] = deal(cell(max(nFam),1)); - [invMZR{1},invOTR{1}] = deal(1); - for nF=2:max(nFam) - n_nF = sum(nFam==nF); - n_nFMZ = sum(nFam(1:nMZF)==nF); + index = ACEcode(hha); + if ( index==ACEcode([0 0 1]) ) + % Parameter A is zero, test-statistic is zero + T = 0; + else + % Parameter A is non-zero... compute log-likelihood... + det_Va = 1; + [invMZR,invOTR] = deal(cell(max(nFam),1)); + [invMZR{1},invOTR{1}] = deal(1); + + for nF=2:max(nFam) + n_nF = sum(nFam==nF); + n_nFMZ = sum(nFam(1:nMZF)==nF); + + %%% AE model + a = hha(1)/2+hha(3); + b = hha(1)/2; + Ra = eye(nF)*a + ones(nF)*b; + det_Va = det_Va*det(Ra)^(n_nF-n_nFMZ); + invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + + Ra(1,2) = hha(1); + Ra(2,1) = Ra(1,2); + det_Va = det_Va*det(Ra)^n_nFMZ; + invMZR{nF} = sparse(pinv(Ra)); + end + + % AE model + Ua = blkdiag(invMZR{nFam(1:nMZF)}, invOTR{nFam(nMZF+1:end)}); + Za = pinv(X'*Ua*X); + Pa = Ua-Ua*X*Za*X'*Ua; + rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; - %%% AE model - a = hha(1)/2+hha(3); - b = hha(1)/2; - Ra = eye(nF)*a + ones(nF)*b; - det_Va = det_Va*det(Ra)^(n_nF-n_nFMZ); - invOTR{nF} = (eye(nF) - ones(nF)*b/(a+b*nF))/a; + % E model + rl0 = -(iY'*P0*iY/sig2+log(det(X'*X)))/2; - Ra(1,2) = hha(1); - Ra(2,1) = Ra(1,2); - det_Va = det_Va*det(Ra)^n_nFMZ; - invMZR{nF} = sparse(pinv(Ra)); + % Compute LRT statistic + T = max(0,-2*(rl0-rla)); end - % AE model - Ua = blkdiag(invMZR{nFam(1:nMZF)}, invOTR{nFam(nMZF+1:end)}); - Za = pinv(X'*Ua*X); - Pa = Ua-Ua*X*Za*X'*Ua; - rla = -(iY'*Pa*iY/sig2+log(det_Va)+log(det(X'*Ua*X)))/2; + Stats(i) = T; - % E model - rl0 = -(iY'*P0*iY/sig2+log(det(X'*X)))/2; - - % Compute LRT statistic - T = max(0,-2*(rl0-rla)); end - - Stats(i) = T; - end end @@ -202,6 +220,10 @@ end +if BadElm>0 + fprintf('WARNING: Voxels/elements found with zero environmental variance (%d). Check for serverely sparse discrete data.\n',BadElm); +end + [mK,mM,mT,uCount] = deal(0); diff --git a/apace/ACEfit_Results.m b/apace/ACEfit_Results.m index 832a81d..458631e 100644 --- a/apace/ACEfit_Results.m +++ b/apace/ACEfit_Results.m @@ -1,4 +1,17 @@ -function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) +function ACEfit_Results(ACEfit_Par,varargin) +% FORMAT ACEfit_Results(ACEfit_Par[,FWEalpha,FDRalpha]) + +if nargin>=2 + FWEalpha = varargin{1}; +else + FWEalpha = 0.05; +end +if nargin>=3 + FDRalpha = varargin{2}; +else + FDRalpha = 0.05; +end + % % Plotting Function % @@ -11,7 +24,7 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) %%% (1) MEAN f = f+1; -figure(f); +SetFig(f); [sMEAN,oMEAN] = sort(mean_ACE); ctl_val_MEAN = sMEAN(ctl_val_index); [f1,x1] = hist(sMEAN,100); @@ -21,29 +34,19 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_MEAN = (N-n_MEAN+1)/N; -set(gca, 'xtick', 0) -if ctl_val_MEAN==mean_ACE(N) - set(gca, 'xtick', ctl_val_MEAN) -else - set(gca, 'xtick', sort([ctl_val_MEAN mean_ACE(N)])) -end - title(sprintf('H0 dist of Mean of h^2, P-value=%.3f',p_MEAN)); xlabel('mean of h^2'); -hold on; yLimits = get(gca,'YLim'); line([mean_ACE(N) mean_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_MEAN ctl_val_MEAN],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); -hold off; -set(gcf,'PaperPosition',[0 0 10 6]) -set(gcf,'PaperSize',[10 6]) +text(mean_ACE(N),0.8*yLimits(2),num2str(mean_ACE(N))) +text(ctl_val_MEAN,0.9*yLimits(2),num2str(ctl_val_MEAN)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_mean.pdf')); - %%% (2) wh2 f = f+1; -figure(f); +SetFig(f); [swh2,owh2] = sort(wh2_ACE); ctl_val_wh2 = swh2(ctl_val_index); [f1,x1] = hist(swh2,100); @@ -53,29 +56,20 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_wh2 = (N-n_wh2+1)/N; -set(gca, 'xtick', 0) -if ctl_val_wh2==wh2_ACE(N) - set(gca, 'xtick', ctl_val_wh2) -else - set(gca, 'xtick', sort([ctl_val_wh2 wh2_ACE(N)])) -end - title(sprintf('H0 dist of Weighted Mean of h^2, P-value=%.3f',p_wh2)); xlabel('weighted mean of h^2'); -hold on; yLimits = get(gca,'YLim'); line([wh2_ACE(N) wh2_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_wh2 ctl_val_wh2],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); -hold off; -set(gcf,'PaperPosition',[0 0 10 6]) -set(gcf,'PaperSize',[10 6]) +text(wh2_ACE(N),0.8*yLimits(2),num2str(wh2_ACE(N))) +text(ctl_val_wh2,0.9*yLimits(2),num2str(ctl_val_wh2)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_wh2.pdf')); %%% (3) MEDIAN f = f+1; -figure(f); +SetFig(f); [sMEDIAN,oMEDIAN] = sort(med_ACE); ctl_val_MEDIAN = sMEDIAN(ctl_val_index); [f1,x1] = hist(sMEDIAN,100); @@ -85,29 +79,20 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_MEDIAN = (N-n_MEDIAN+1)/N; -set(gca, 'xtick', 0) -if ctl_val_MEDIAN==med_ACE(N) - set(gca, 'xtick', ctl_val_MEDIAN) -else - set(gca, 'xtick', sort([ctl_val_MEDIAN med_ACE(N)])) -end - title(sprintf('H0 dist of Median (Q2) of h^2, P-value=%.3f',p_MEDIAN)); xlabel('Q2 of h^2'); -hold on; yLimits = get(gca,'YLim'); line([med_ACE(N) med_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_MEDIAN ctl_val_MEDIAN],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); -hold off; -set(gcf,'PaperPosition',[0 0 10 6]) -set(gcf,'PaperSize',[10 6]) +text(med_ACE(N),0.8*yLimits(2),num2str(med_ACE(N))) +text(ctl_val_MEDIAN,0.9*yLimits(2),num2str(ctl_val_MEDIAN)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_median.pdf')); %%% (4) Q3 f = f+1; -figure(f); +SetFig(f); [sQ3,oQ3] = sort(q3_ACE); ctl_val_Q3 = sQ3(ctl_val_index); [f1,x1] = hist(sQ3,100); @@ -117,23 +102,14 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_Q3 = (N-n_Q3+1)/N; -set(gca, 'xtick', 0) -if ctl_val_Q3==q3_ACE(N) - set(gca, 'xtick', ctl_val_Q3) -else - set(gca, 'xtick', sort([ctl_val_Q3 q3_ACE(N)])) -end - title(sprintf('H0 dist of Third Quartile (Q3) of h^2, P-value=%.3f',p_Q3)); xlabel('Q3 of h^2'); -hold on; yLimits = get(gca,'YLim'); line([q3_ACE(N) q3_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_Q3 ctl_val_Q3],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); -hold off; -set(gcf,'PaperPosition',[0 0 10 6]) -set(gcf,'PaperSize',[10 6]) +text(q3_ACE(N),0.8*yLimits(2),num2str(q3_ACE(N))) +text(ctl_val_Q3,0.9*yLimits(2),num2str(ctl_val_Q3)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_q3.pdf')); @@ -145,7 +121,7 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) else f = f+1; - figure(f); + SetFig(f); [smGTmedian,omGTmedian] = sort(mGmed_ACE); ctl_val_mGTmedian = smGTmedian(ctl_val_index); [f1,x1] = hist(smGTmedian,100); @@ -155,25 +131,16 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_mGTmedian = (N-n_mGTmedian+1)/N; - set(gca, 'xtick', 0) - if ctl_val_mGTmedian==mGmed_ACE(N) - set(gca, 'xtick', ctl_val_mGTmedian) - else - set(gca, 'xtick', sort([ctl_val_mGTmedian mGmed_ACE(N)])) - end - title(sprintf('H0 dist of Mean of h^2 > Q2(h^2), P-value=%.3f',p_mGTmedian)); xlabel('mean of h^2 > Q2(h^2)'); - hold on; yLimits = get(gca,'YLim'); line([mGmed_ACE(N) mGmed_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_mGTmedian ctl_val_mGTmedian],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(mGmed_ACE(N),0.8*yLimits(2),num2str(mGmed_ACE(N))) + text(ctl_val_mGTmedian,0.9*yLimits(2),num2str(ctl_val_mGTmedian)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_mGTmedian.pdf')); - + end @@ -185,7 +152,7 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) else f = f+1; - figure(f); + SetFig(f); [smGTq3,omGTq3] = sort(mGq3_ACE); ctl_val_mGTq3 = smGTq3(ctl_val_index); [f1,x1] = hist(smGTq3,100); @@ -195,23 +162,14 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_mGTq3 = (N-n_mGTq3+1)/N; - set(gca, 'xtick', 0) - if ctl_val_mGTq3==mGq3_ACE(N) - set(gca, 'xtick', ctl_val_mGTq3) - else - set(gca, 'xtick', sort([ctl_val_mGTq3 mGq3_ACE(N)])) - end - title(sprintf('H0 dist of Mean of h^2 > Q3(h^2), P-value=%.3f',p_mGTq3)); xlabel('mean of h^2 > Q3(h^2)'); - hold on; yLimits = get(gca,'YLim'); line([mGq3_ACE(N) mGq3_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_mGTq3 ctl_val_mGTq3],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(mGq3_ACE(N),0.8*yLimits(2),num2str(mGq3_ACE(N))) + text(ctl_val_mGTq3,0.9*yLimits(2),num2str(ctl_val_mGTq3)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_mGTq3.pdf')); end @@ -226,7 +184,7 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) %%% (7) Maximum Test Statistic T f = f+1; - figure(f); + SetFig(f); [sT,oT] = sort(max_T_ACE); ctl_val_T = sT(ctl_val_index); [f3,x3] = hist(sT,100); @@ -236,24 +194,15 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_T = (N-n_T+1)/N; - set(gca, 'xtick', 0) - if ctl_val_T==max_T_ACE(N) - set(gca, 'xtick', ctl_val_T) - else - set(gca, 'xtick', sort([ctl_val_T max_T_ACE(N)])) - end - title(sprintf('H0 dist of Maximum LRT Statistic (T), FWE P-value=%.3f',p_T)); xlabel('T'); ylabel('f(T)'); - hold on; yLimits = get(gca,'YLim'); line([max_T_ACE(N) max_T_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_T ctl_val_T],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(max_T_ACE(N),0.8*yLimits(2),num2str(max_T_ACE(N))) + text(ctl_val_T,0.9*yLimits(2),num2str(ctl_val_T)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_test_statistic.pdf')); fprintf('The %.2f level critical threshold for the maximum LRT statistic is %8.3f. \n \n', FWEalpha, ctl_val_T); @@ -287,7 +236,7 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) %%% (8) FDR plot f = f+1; - figure(f); + SetFig(f); plot((1:nMsk)'/nMsk,sort(unPval_ACE(ACEfit_Par.I_data')),(1:nMsk)'/nMsk,sort(fdrPval)); legend({'Uncorr P','FDR-corr P'},'Location','NorthWest') xlabel('Expected Uncorrected Ordered P'); @@ -296,8 +245,6 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) axis equal;axis tight;axis([0 1 0 1]); abline(0,1,'LineStyle','-','color',[.5 .5 .5]); abline(0,FDRalpha,'LineStyle',':','color','red') - set(gcf,'PaperPosition',[0 0 6 6]) - set(gcf,'PaperSize',[6 6]) print('-dpdf',fullfile(ACEfit_Par.ResDir,'h2_FDRplot.pdf')); FDR_Pval = ones(Dim); @@ -329,19 +276,26 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) end fprintf('Minimum element-wise P_FWE is %.4f. \n \n', 10.^(-max(corrPval_ACE))); - + % % Write out the output images % - WriteData(unPval_ACE, Vs, 'ACE_A_LRT_vox_P', ACEfit_Par.ResDir); - WriteData(corrPval_ACE, Vs, 'ACE_A_LRT_vox_FWEP', ACEfit_Par.ResDir); - WriteData(FDR_Pval, Vs, 'ACE_A_LRT_vox_FDRP', ACEfit_Par.ResDir); + switch upper(ACEfit_Par.Model) + case 'ACE' + WriteData(unPval_ACE, Vs, 'ACE_A_LRT_vox_P', ACEfit_Par.ResDir); + WriteData(corrPval_ACE, Vs, 'ACE_A_LRT_vox_FWEP', ACEfit_Par.ResDir); + WriteData(FDR_Pval, Vs, 'ACE_A_LRT_vox_FDRP', ACEfit_Par.ResDir); + case 'AE' + WriteData(unPval_ACE, Vs, 'AE_A_LRT_vox_P', ACEfit_Par.ResDir); + WriteData(corrPval_ACE, Vs, 'AE_A_LRT_vox_FWEP', ACEfit_Par.ResDir); + WriteData(FDR_Pval, Vs, 'AE_A_LRT_vox_FDRP', ACEfit_Par.ResDir); + end if ACEfit_Par.Vs.ClustInf %%% (9) Maximum Suprathreshold Cluster Size K f = f+1; - figure(f); + SetFig(f); [sK,oK] = sort(max_K_ACE); ctl_val_K = sK(ctl_val_index); [f1,x1] = hist(sK,100); @@ -351,30 +305,21 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_K = (N-n_K+1)/N; - set(gca, 'xtick', 0) - if ctl_val_K==max_K_ACE(N) - set(gca, 'xtick', ctl_val_K) - else - set(gca, 'xtick', sort([ctl_val_K max_K_ACE(N)])) - end - title(sprintf('H0 dist of Maximum Suprathreshold Cluster Size (K), FWE P-value=%.3f',p_K)); xlabel('K'); ylabel('f(K)'); - hold on; yLimits = get(gca,'YLim'); line([max_K_ACE(N) max_K_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_K ctl_val_K],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(max_K_ACE(N),0.8*yLimits(2),num2str(max_K_ACE(N))) + text(ctl_val_K,0.9*yLimits(2),num2str(ctl_val_K)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_cluster_size.pdf')); %%% (10) Maximum Suprathreshold Cluster Mass M f = f+1; - figure(f); + SetFig(f); [sM,oM] = sort(max_M_ACE); ctl_val_M = sM(ctl_val_index); [f2,x2] = hist(sM,100); @@ -384,24 +329,15 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % Calculate the permutation-based p-value p_M = (N-n_M+1)/N; - set(gca, 'xtick', 0) - if ctl_val_M==max_M_ACE(N) - set(gca, 'xtick', ctl_val_M) - else - set(gca, 'xtick', sort([ctl_val_M max_M_ACE(N)])) - end - title(sprintf('H0 dist of Maximum Suprathreshold Cluster Mass (M), FWE P-value = %.3f',p_M)); xlabel('M'); ylabel('f(M)'); - hold on; yLimits = get(gca,'YLim'); line([max_M_ACE(N) max_M_ACE(N)],[0 yLimits(2)],'Marker','.','Color','green'); % Plot the critical threshold (red) line([ctl_val_M ctl_val_M],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(max_M_ACE(N),0.8*yLimits(2),num2str(max_M_ACE(N))) + text(ctl_val_M,0.9*yLimits(2),num2str(ctl_val_M)) print('-dpdf',fullfile(ACEfit_Par.ResDir,'H0dist_cluster_mass.pdf')); fprintf('The %.2f level critical threshold for the maximum suprathreshold cluster size is %d. \n', FWEalpha, ctl_val_K); @@ -441,10 +377,18 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) % % Write out the output images % - WriteData(ClSz, Vs, 'ACE_A_LRT_clus', ACEfit_Par.ResDir); - WriteData(ClMass, Vs, 'ACE_A_LRT_mass', ACEfit_Par.ResDir); - WriteData(corrPClSz, Vs, 'ACE_A_LRT_clus_FWEP', ACEfit_Par.ResDir); - WriteData(corrPClMass, Vs, 'ACE_A_LRT_mass_FWEP', ACEfit_Par.ResDir); + switch upper(ACEfit_Par.Model) + case 'ACE' + WriteData(ClSz, Vs, 'ACE_A_LRT_clus', ACEfit_Par.ResDir); + WriteData(ClMass, Vs, 'ACE_A_LRT_mass', ACEfit_Par.ResDir); + WriteData(corrPClSz, Vs, 'ACE_A_LRT_clus_FWEP', ACEfit_Par.ResDir); + WriteData(corrPClMass, Vs, 'ACE_A_LRT_mass_FWEP', ACEfit_Par.ResDir); + case 'AE' + WriteData(ClSz, Vs, 'AE_A_LRT_clus', ACEfit_Par.ResDir); + WriteData(ClMass, Vs, 'AE_A_LRT_mass', ACEfit_Par.ResDir); + WriteData(corrPClSz, Vs, 'AE_A_LRT_clus_FWEP', ACEfit_Par.ResDir); + WriteData(corrPClMass, Vs, 'AE_A_LRT_mass_FWEP', ACEfit_Par.ResDir); + end save(fullfile(ACEfit_Par.ResDir,'Pvals_Max_h2'),'p_T','p_K','p_M'); @@ -456,4 +400,12 @@ function ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha) end -return \ No newline at end of file +return + +function SetFig(f) + +figure(f) +set(gcf,'PaperPosition',[0 0 10 6]) +set(gcf,'PaperSize',[10 6]) + +return diff --git a/apace/APACEsummary.m b/apace/APACEsummary.m index ba7e864..f156c0c 100644 --- a/apace/APACEsummary.m +++ b/apace/APACEsummary.m @@ -23,7 +23,7 @@ Names = {'Mean','VarWtMean','AgHe','Median Q2','Q3','Mean(>Q2)','Mean(>Q3)'}; % Names = {'Mean','VarWtMean','AgHeNorm','AgHeNoNorm','Median Q2','Q3','Mean(>Q2)','Mean(>Q3)'}; -if exist(fullfile(ACEfit_Par.ResDir,'Ests_AgHe'),'file') +if exist(fullfile(ACEfit_Par.ResDir,'Ests_AgHe.mat'),'file') E_s = load(fullfile(ACEfit_Par.ResDir,'Ests_AgHe')); % Est_MZDZ Est_DZSib % Ne = load(fullfile(ACEfit_Par.ResDir,'Ests_AgHe_Norm')); % Oe = load(fullfile(ACEfit_Par.ResDir,'Ests_AgHe_NoNorm')); @@ -76,7 +76,7 @@ load(fullfile(ACEfit_Par.ResDir,'Pvals_h2')) % Pvals_h2 - if exist(fullfile(ACEfit_Par.ResDir,'Pvals_AgHe'),'file') + if exist(fullfile(ACEfit_Par.ResDir,'Pvals_AgHe.mat'),'file') P_s = load(fullfile(ACEfit_Par.ResDir,'Pvals_AgHe')); % Pvals_MZDZ Pvals_DZSib % Np = load(fullfile(ACEfit_Par.ResDir,'Pvals_AgHe_Norm')); % Op = load(fullfile(ACEfit_Par.ResDir,'Pvals_AgHe_NoNorm')); @@ -111,7 +111,7 @@ load(fullfile(ACEfit_Par.ResDir,'Boot_CIs')) % alpha CIs_h2 CIs_c2 CIs_e2 - if exist(fullfile(ACEfit_Par.ResDir,'CIs_AgHe'),'file') + if exist(fullfile(ACEfit_Par.ResDir,'CIs_AgHe.mat'),'file') CI_s = load(fullfile(ACEfit_Par.ResDir,'CIs_AgHe')); % CI_MZDZ CI_DZSib % Nc = load(fullfile(ACEfit_Par.ResDir,'CIs_AgHe_Norm')); % Oc = load(fullfile(ACEfit_Par.ResDir,'CIs_AgHe_NoNorm')); @@ -177,7 +177,7 @@ CIs = [ CIs; {[NaN NaN],[NaN NaN],[NaN NaN]} ]; - if ACEfit_Par.nPerm>0 + if ACEfit_Par.nPerm>0 load(fullfile(ACEfit_Par.ResDir,'Pvals_Max_h2')) % FWE voxel-wise p-value p_T (and cluster-based p-value p_K & p_M) Ps = [ Ps; p_T, NaN, NaN ]; diff --git a/apace/AgHe_Method.m b/apace/AgHe_Method.m index bab0fd4..b9401d5 100644 --- a/apace/AgHe_Method.m +++ b/apace/AgHe_Method.m @@ -1,21 +1,32 @@ -function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) +function AgHe_Method(ACEfit_Par,varargin) % % Permutation & bootstrap inferences for AgHe (aka Steve's) method % -% AgHe_Method(ACEfit_Par,Palpha,Balpha) - Save results -% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_Norm') - Save results with the suffix -% +% AgHe_Method(ACEfit_Par) - Save results +% AgHe_Method(ACEfit_Par,'_Norm') - Save results with suffix +% for permutation & bootstrap +% AgHe_Method(ACEfit_Par,'_Norm',Palpha,Balpha) - Save results with the suffix, +% also set Perm & Bootstrap alpha if size(ACEfit_Par.Y,1)==1 error('Aggregate heritability is designed for multiple phenotypes!') end -if nargin<=3 +if nargin<2 ResSuf = ''; else - % Results suffix ResSuf = varargin{1}; end +if nargin<3 || isempty(varargin{2}) + Palpha = 0.05; +else + Palpha = varargin{2}; +end +if nargin<4 || isempty(varargin{3}) + Balpha = 0.05; +else + Balpha = varargin{3}; +end nMZF = ACEfit_Par.nMZF; nDZF = ACEfit_Par.nDZF; @@ -30,6 +41,10 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % Res = ACEfit_Par.Y' - ACEfit_Par.X*pinv(ACEfit_Par.X)*ACEfit_Par.Y'; +if ( ~isfield(ACEfit_Par,'AggNlz') || isempty(ACEfit_Par.AggNlz) ) + ACEfit_Par.AggNlz = 0; +end + switch(ACEfit_Par.AggNlz) case 1 for j=1:size(Res,2) @@ -149,7 +164,6 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) Est_DZSib = mean(rDZo) - mean(rSibo); - f = get(0,'Children'); if isempty(f) f = 1; @@ -158,7 +172,7 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % f = max(f)+1; end -figure(f); +SetFig(f); % AgHe box plots Z = [rMZo; rDZo; rSibo ]; g = [zeros(length(rMZo),1); ones(length(rDZo),1); 2*ones(length(rSibo),1)]; @@ -166,16 +180,14 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) title(sprintf('Boxplots of Correlations of Aggregate Heritabilities')); set(gca, 'XTick', [1; 2; 3]); set(gca, 'XTickLabel', {'MZ'; 'DZ'; 'SIB'}); -set(gcf,'PaperPosition',[0 0 10 6]) -set(gcf,'PaperSize',[10 6]) print('-dpdf',fullfile(ACEfit_Par.ResDir,'AgHe_CorrPlots.pdf')); save(fullfile(ACEfit_Par.ResDir,['Ests_AgHe' ResSuf]),... - 'Est_MZDZ','Est_DZSib',... - 'rMZo', 'rDZo', 'rSibo'); + 'Est_MZDZ','Est_DZSib',... + 'rMZo', 'rDZo', 'rSibo'); % @@ -236,7 +248,7 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % 1st statistic: 2-sample t-statistic % f = f+1; - figure(f); + SetFig(f); [st_stats_MZ_DZ,~] = sort(t_stats_MZ_DZ); ctl_val_t_stats_MZ_DZ = st_stats_MZ_DZ(ctl_val_index); [f1,x1] = hist(st_stats_MZ_DZ,100); @@ -246,28 +258,19 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % Calculate the permutation-based p-value p_t_stats_MZ_DZ = (N-n_t_stats_MZ_DZ+1)/N; - set(gca, 'xtick', 0) - if ctl_val_t_stats_MZ_DZ==t_stats_MZ_DZ(N) - set(gca, 'xtick', ctl_val_t_stats_MZ_DZ) - else - set(gca, 'xtick', sort([ctl_val_t_stats_MZ_DZ t_stats_MZ_DZ(N)])) - end - title(sprintf('H0 dist of two-sample t-statistic for rMZ & rDZ, P-value=%.3f',p_t_stats_MZ_DZ)); xlabel('t-statistic for rMZ vs. rDZ'); ylabel('frequency'); - hold on; yLimits = get(gca,'YLim'); line([t_stats_MZ_DZ(N) t_stats_MZ_DZ(N)],[0 yLimits(2)],'Marker','.','Color','green'); line([ctl_val_t_stats_MZ_DZ ctl_val_t_stats_MZ_DZ],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(t_stats_MZ_DZ(N),0.8*yLimits(2),num2str(t_stats_MZ_DZ(N))) + text(ctl_val_t_stats_MZ_DZ,0.9*yLimits(2),num2str(ctl_val_t_stats_MZ_DZ)) print('-dpdf',fullfile(ACEfit_Par.ResDir,['H0dist_Tstat_rMZ_rDZ' ResSuf '.pdf'])); f = f+1; - figure(f); + SetFig(f); [st_stats_DZ_Sib,~] = sort(t_stats_DZ_Sib); ctl_val_t_stats_DZ_Sib = st_stats_DZ_Sib(ctl_val_index); [f1,x1] = hist(st_stats_DZ_Sib,100); @@ -277,23 +280,14 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % Calculate the permutation-based p-value p_t_stats_DZ_Sib = (N-n_t_stats_DZ_Sib+1)/N; - set(gca, 'xtick', 0) - if ctl_val_t_stats_DZ_Sib==t_stats_DZ_Sib(N) - set(gca, 'xtick', ctl_val_t_stats_DZ_Sib) - else - set(gca, 'xtick', sort([ctl_val_t_stats_DZ_Sib t_stats_DZ_Sib(N)])) - end - title(sprintf('H0 dist of two-sample t-statistic for rDZ & rSib, P-value=%.3f',p_t_stats_DZ_Sib)); xlabel('t-statistic for rDZ vs. rSib'); ylabel('frequency'); - hold on; yLimits = get(gca,'YLim'); line([t_stats_DZ_Sib(N) t_stats_DZ_Sib(N)],[0 yLimits(2)],'Marker','.','Color','green'); line([ctl_val_t_stats_DZ_Sib ctl_val_t_stats_DZ_Sib],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(t_stats_DZ_Sib(N),0.8*yLimits(2),num2str(t_stats_DZ_Sib(N))) + text(ctl_val_t_stats_DZ_Sib,0.9*yLimits(2),num2str(ctl_val_t_stats_DZ_Sib)) print('-dpdf',fullfile(ACEfit_Par.ResDir,['H0dist_Tstat_rDZ_rSib' ResSuf '.pdf'])); @@ -301,7 +295,7 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % 2nd statistic: mean difference % f = f+1; - figure(f); + SetFig(f); [sCorrDiff_MZ_DZ,~] = sort(CorrDiff_MZ_DZ); ctl_val_CorrDiff_MZ_DZ = sCorrDiff_MZ_DZ(ctl_val_index); [f2,x2] = hist(sCorrDiff_MZ_DZ,100); @@ -311,28 +305,19 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % Calculate the permutation-based p-value p_CorrDiff_MZ_DZ = (N-n_CorrDiff_MZ_DZ+1)/N; - set(gca, 'xtick', 0) - if ctl_val_CorrDiff_MZ_DZ==CorrDiff_MZ_DZ(N) - set(gca, 'xtick', ctl_val_CorrDiff_MZ_DZ) - else - set(gca, 'xtick', sort([ctl_val_CorrDiff_MZ_DZ CorrDiff_MZ_DZ(N)])) - end - title(sprintf('H0 dist of mean difference between rMZ & rDZ, P-value=%.3f',p_CorrDiff_MZ_DZ)); xlabel('E(rMZ) - E(rDZ)'); ylabel('frequency'); - hold on; yLimits = get(gca,'YLim'); line([CorrDiff_MZ_DZ(N) CorrDiff_MZ_DZ(N)],[0 yLimits(2)],'Marker','.','Color','green'); line([ctl_val_CorrDiff_MZ_DZ ctl_val_CorrDiff_MZ_DZ],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(CorrDiff_MZ_DZ(N),0.8*yLimits(2),num2str(CorrDiff_MZ_DZ(N))) + text(ctl_val_CorrDiff_MZ_DZ,0.9*yLimits(2),num2str(ctl_val_CorrDiff_MZ_DZ)) print('-dpdf',fullfile(ACEfit_Par.ResDir,['H0dist_Diff_rMZ_rDZ' ResSuf '.pdf'])); f = f+1; - figure(f); + SetFig(f); [sCorrDiff_DZ_Sib,~] = sort(CorrDiff_DZ_Sib); ctl_val_CorrDiff_DZ_Sib = sCorrDiff_DZ_Sib(ctl_val_index); [f2,x2] = hist(sCorrDiff_DZ_Sib,100); @@ -342,23 +327,14 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) % Calculate the permutation-based p-value p_CorrDiff_DZ_Sib = (N-n_CorrDiff_DZ_Sib+1)/N; - set(gca, 'xtick', 0) - if ctl_val_CorrDiff_DZ_Sib==CorrDiff_DZ_Sib(N) - set(gca, 'xtick', ctl_val_CorrDiff_DZ_Sib) - else - set(gca, 'xtick', sort([ctl_val_CorrDiff_DZ_Sib CorrDiff_DZ_Sib(N)])) - end - title(sprintf('H0 dist of mean difference between rDZ & rSib, P-value=%.3f',p_CorrDiff_DZ_Sib)); xlabel('E[rDZ-rSib]'); ylabel('frequency'); - hold on; yLimits = get(gca,'YLim'); line([CorrDiff_DZ_Sib(N) CorrDiff_DZ_Sib(N)],[0 yLimits(2)],'Marker','.','Color','green'); line([ctl_val_CorrDiff_DZ_Sib ctl_val_CorrDiff_DZ_Sib],[0 yLimits(2)],'Marker','.','LineStyle','-.','Color','red'); - hold off; - set(gcf,'PaperPosition',[0 0 10 6]) - set(gcf,'PaperSize',[10 6]) + text(CorrDiff_DZ_Sib(N),0.8*yLimits(2),num2str(CorrDiff_DZ_Sib(N))) + text(ctl_val_CorrDiff_DZ_Sib,0.9*yLimits(2),num2str(ctl_val_CorrDiff_DZ_Sib)) print('-dpdf',fullfile(ACEfit_Par.ResDir,['H0dist_Diff_rDZ_rSib' ResSuf '.pdf'])); @@ -423,5 +399,10 @@ function AgHe_Method(ACEfit_Par,Palpha,Balpha,varargin) return +function SetFig(f) +figure(f) +set(gcf,'PaperPosition',[0 0 10 6]) +set(gcf,'PaperSize',[10 6]) +return diff --git a/apace/Boot_CIs.m b/apace/Boot_CIs.m index a0003fb..64b8803 100644 --- a/apace/Boot_CIs.m +++ b/apace/Boot_CIs.m @@ -1,8 +1,15 @@ -function Boot_CIs(ACEfit_Par,alpha) +function Boot_CIs(ACEfit_Par,varargin) +% FORMAT Boot_CIs(ACEfit_Par[,Balpha]) % % Generate Bootstrap CIs. % +if nargin>=2 + Balpha = varargin{1}; +else + Balpha = 0.05; +end + load(fullfile(ACEfit_Par.ResDir,'ACEfit_Boot')); switch upper(ACEfit_Par.Model) @@ -14,32 +21,32 @@ function Boot_CIs(ACEfit_Par,alpha) % % (1) mean, wh2, median, q3, mean(h2>median), mean(h2>q3) % - [CIs_h2(1,1), CIs_h2(1,2)] = Bootstrap_CI(meanh2_ACE, alpha); - [CIs_h2(2,1), CIs_h2(2,2)] = Bootstrap_CI(wh2_ACE, alpha); - [CIs_h2(3,1), CIs_h2(3,2)] = Bootstrap_CI(medh2_ACE, alpha); - [CIs_h2(4,1), CIs_h2(4,2)] = Bootstrap_CI(q3h2_ACE, alpha); - [CIs_h2(5,1), CIs_h2(5,2)] = Bootstrap_CI(mGmedh2_ACE, alpha); - [CIs_h2(6,1), CIs_h2(6,2)] = Bootstrap_CI(mGq3h2_ACE, alpha); + [CIs_h2(1,1), CIs_h2(1,2)] = Bootstrap_CI(meanh2_ACE, Balpha); + [CIs_h2(2,1), CIs_h2(2,2)] = Bootstrap_CI(wh2_ACE, Balpha); + [CIs_h2(3,1), CIs_h2(3,2)] = Bootstrap_CI(medh2_ACE, Balpha); + [CIs_h2(4,1), CIs_h2(4,2)] = Bootstrap_CI(q3h2_ACE, Balpha); + [CIs_h2(5,1), CIs_h2(5,2)] = Bootstrap_CI(mGmedh2_ACE, Balpha); + [CIs_h2(6,1), CIs_h2(6,2)] = Bootstrap_CI(mGq3h2_ACE, Balpha); % % (2) mean, wc2, median, q3, mean(c2>median), mean(c2>q3) % - [CIs_c2(1,1), CIs_c2(1,2)] = Bootstrap_CI(meanc2_ACE, alpha); - [CIs_c2(2,1), CIs_c2(2,2)] = Bootstrap_CI(wc2_ACE, alpha); - [CIs_c2(3,1), CIs_c2(3,2)] = Bootstrap_CI(medc2_ACE, alpha); - [CIs_c2(4,1), CIs_c2(4,2)] = Bootstrap_CI(q3c2_ACE, alpha); - [CIs_c2(5,1), CIs_c2(5,2)] = Bootstrap_CI(mGmedc2_ACE, alpha); - [CIs_c2(6,1), CIs_c2(6,2)] = Bootstrap_CI(mGq3c2_ACE, alpha); + [CIs_c2(1,1), CIs_c2(1,2)] = Bootstrap_CI(meanc2_ACE, Balpha); + [CIs_c2(2,1), CIs_c2(2,2)] = Bootstrap_CI(wc2_ACE, Balpha); + [CIs_c2(3,1), CIs_c2(3,2)] = Bootstrap_CI(medc2_ACE, Balpha); + [CIs_c2(4,1), CIs_c2(4,2)] = Bootstrap_CI(q3c2_ACE, Balpha); + [CIs_c2(5,1), CIs_c2(5,2)] = Bootstrap_CI(mGmedc2_ACE, Balpha); + [CIs_c2(6,1), CIs_c2(6,2)] = Bootstrap_CI(mGq3c2_ACE, Balpha); % % (3) mean, we2, median, q3, mean(e2>median), mean(e2>q3) % - [CIs_e2(1,1), CIs_e2(1,2)] = Bootstrap_CI(meane2_ACE, alpha); - [CIs_e2(2,1), CIs_e2(2,2)] = Bootstrap_CI(we2_ACE, alpha); - [CIs_e2(3,1), CIs_e2(3,2)] = Bootstrap_CI(mede2_ACE, alpha); - [CIs_e2(4,1), CIs_e2(4,2)] = Bootstrap_CI(q3e2_ACE, alpha); - [CIs_e2(5,1), CIs_e2(5,2)] = Bootstrap_CI(mGmede2_ACE, alpha); - [CIs_e2(6,1), CIs_e2(6,2)] = Bootstrap_CI(mGq3e2_ACE, alpha); + [CIs_e2(1,1), CIs_e2(1,2)] = Bootstrap_CI(meane2_ACE, Balpha); + [CIs_e2(2,1), CIs_e2(2,2)] = Bootstrap_CI(we2_ACE, Balpha); + [CIs_e2(3,1), CIs_e2(3,2)] = Bootstrap_CI(mede2_ACE, Balpha); + [CIs_e2(4,1), CIs_e2(4,2)] = Bootstrap_CI(q3e2_ACE, Balpha); + [CIs_e2(5,1), CIs_e2(5,2)] = Bootstrap_CI(mGmede2_ACE, Balpha); + [CIs_e2(6,1), CIs_e2(6,2)] = Bootstrap_CI(mGq3e2_ACE, Balpha); CIs_h2 = max(0,CIs_h2); CIs_h2 = min(1,CIs_h2); @@ -50,28 +57,28 @@ function Boot_CIs(ACEfit_Par,alpha) CIs_e2 = max(0,CIs_e2); CIs_e2 = min(1,CIs_e2); - fprintf('The %.1f%% confidence interval for mean of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(1,:)); - fprintf('The %.1f%% confidence interval for weighted mean of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(2,:)); - fprintf('The %.1f%% confidence interval for median (Q2) of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(3,:)); - fprintf('The %.1f%% confidence interval for the third quantile (Q3) of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(4,:)); - fprintf('The %.1f%% confidence interval for mean of h2>Q2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(5,:)); - fprintf('The %.1f%% confidence interval for mean of h2>Q3 is (%.2f, %.2f). \n \n', 100*(1-alpha), CIs_h2(6,:)); - - fprintf('The %.1f%% confidence interval for mean of c2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_c2(1,:)); - fprintf('The %.1f%% confidence interval for weighted mean of c2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_c2(2,:)); - fprintf('The %.1f%% confidence interval for median (Q2) of c2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_c2(3,:)); - fprintf('The %.1f%% confidence interval for the third quantile (Q3) of c2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_c2(4,:)); - fprintf('The %.1f%% confidence interval for mean of c2>Q2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_c2(5,:)); - fprintf('The %.1f%% confidence interval for mean of c2>Q3 is (%.2f, %.2f). \n \n', 100*(1-alpha), CIs_c2(6,:)); - - fprintf('The %.1f%% confidence interval for mean of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(1,:)); - fprintf('The %.1f%% confidence interval for weighted mean of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(2,:)); - fprintf('The %.1f%% confidence interval for median (Q2) of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(3,:)); - fprintf('The %.1f%% confidence interval for the third quantile (Q3) of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(4,:)); - fprintf('The %.1f%% confidence interval for mean of e2>Q2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(5,:)); - fprintf('The %.1f%% confidence interval for mean of e2>Q3 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(6,:)); - - save(fullfile(ACEfit_Par.ResDir,'Boot_CIs'),'alpha','CIs_h2','CIs_c2','CIs_e2'); + fprintf('The %.1f%% confidence interval for mean of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(1,:)); + fprintf('The %.1f%% confidence interval for weighted mean of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(2,:)); + fprintf('The %.1f%% confidence interval for median (Q2) of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(3,:)); + fprintf('The %.1f%% confidence interval for the third quantile (Q3) of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(4,:)); + fprintf('The %.1f%% confidence interval for mean of h2>Q2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(5,:)); + fprintf('The %.1f%% confidence interval for mean of h2>Q3 is (%.2f, %.2f). \n \n', 100*(1-Balpha), CIs_h2(6,:)); + + fprintf('The %.1f%% confidence interval for mean of c2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_c2(1,:)); + fprintf('The %.1f%% confidence interval for weighted mean of c2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_c2(2,:)); + fprintf('The %.1f%% confidence interval for median (Q2) of c2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_c2(3,:)); + fprintf('The %.1f%% confidence interval for the third quantile (Q3) of c2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_c2(4,:)); + fprintf('The %.1f%% confidence interval for mean of c2>Q2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_c2(5,:)); + fprintf('The %.1f%% confidence interval for mean of c2>Q3 is (%.2f, %.2f). \n \n', 100*(1-Balpha), CIs_c2(6,:)); + + fprintf('The %.1f%% confidence interval for mean of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(1,:)); + fprintf('The %.1f%% confidence interval for weighted mean of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(2,:)); + fprintf('The %.1f%% confidence interval for median (Q2) of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(3,:)); + fprintf('The %.1f%% confidence interval for the third quantile (Q3) of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(4,:)); + fprintf('The %.1f%% confidence interval for mean of e2>Q2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(5,:)); + fprintf('The %.1f%% confidence interval for mean of e2>Q3 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(6,:)); + + save(fullfile(ACEfit_Par.ResDir,'Boot_CIs'),'Balpha','CIs_h2','CIs_c2','CIs_e2'); case 'AE' @@ -80,22 +87,22 @@ function Boot_CIs(ACEfit_Par,alpha) % % (1) meanh2, wh2, median, q3, mean(h2>median), mean(h2>q3) % - [CIs_h2(1,1), CIs_h2(1,2)] = Bootstrap_CI(meanh2_ACE, alpha); - [CIs_h2(2,1), CIs_h2(2,2)] = Bootstrap_CI(wh2_ACE, alpha); - [CIs_h2(3,1), CIs_h2(3,2)] = Bootstrap_CI(medh2_ACE, alpha); - [CIs_h2(4,1), CIs_h2(4,2)] = Bootstrap_CI(q3h2_ACE, alpha); - [CIs_h2(5,1), CIs_h2(5,2)] = Bootstrap_CI(mGmedh2_ACE, alpha); - [CIs_h2(6,1), CIs_h2(6,2)] = Bootstrap_CI(mGq3h2_ACE, alpha); + [CIs_h2(1,1), CIs_h2(1,2)] = Bootstrap_CI(meanh2_ACE, Balpha); + [CIs_h2(2,1), CIs_h2(2,2)] = Bootstrap_CI(wh2_ACE, Balpha); + [CIs_h2(3,1), CIs_h2(3,2)] = Bootstrap_CI(medh2_ACE, Balpha); + [CIs_h2(4,1), CIs_h2(4,2)] = Bootstrap_CI(q3h2_ACE, Balpha); + [CIs_h2(5,1), CIs_h2(5,2)] = Bootstrap_CI(mGmedh2_ACE, Balpha); + [CIs_h2(6,1), CIs_h2(6,2)] = Bootstrap_CI(mGq3h2_ACE, Balpha); % % (2) meane2, we2, median, q3, mean(e2>median), mean(e2>q3) % - [CIs_e2(1,1), CIs_e2(1,2)] = Bootstrap_CI(meane2_ACE, alpha); - [CIs_e2(2,1), CIs_e2(2,2)] = Bootstrap_CI(we2_ACE, alpha); - [CIs_e2(3,1), CIs_e2(3,2)] = Bootstrap_CI(mede2_ACE, alpha); - [CIs_e2(4,1), CIs_e2(4,2)] = Bootstrap_CI(q3e2_ACE, alpha); - [CIs_e2(5,1), CIs_e2(5,2)] = Bootstrap_CI(mGmede2_ACE, alpha); - [CIs_e2(6,1), CIs_e2(6,2)] = Bootstrap_CI(mGq3e2_ACE, alpha); + [CIs_e2(1,1), CIs_e2(1,2)] = Bootstrap_CI(meane2_ACE, Balpha); + [CIs_e2(2,1), CIs_e2(2,2)] = Bootstrap_CI(we2_ACE, Balpha); + [CIs_e2(3,1), CIs_e2(3,2)] = Bootstrap_CI(mede2_ACE, Balpha); + [CIs_e2(4,1), CIs_e2(4,2)] = Bootstrap_CI(q3e2_ACE, Balpha); + [CIs_e2(5,1), CIs_e2(5,2)] = Bootstrap_CI(mGmede2_ACE, Balpha); + [CIs_e2(6,1), CIs_e2(6,2)] = Bootstrap_CI(mGq3e2_ACE, Balpha); CIs_h2 = max(0,CIs_h2); CIs_h2 = min(1,CIs_h2); @@ -103,21 +110,21 @@ function Boot_CIs(ACEfit_Par,alpha) CIs_e2 = max(0,CIs_e2); CIs_e2 = min(1,CIs_e2); - fprintf('The %.1f%% confidence interval for mean of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(1,:)); - fprintf('The %.1f%% confidence interval for weighted mean of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(2,:)); - fprintf('The %.1f%% confidence interval for median (Q2) of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(3,:)); - fprintf('The %.1f%% confidence interval for the third quantile (Q3) of h2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(4,:)); - fprintf('The %.1f%% confidence interval for mean of h2>Q2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_h2(5,:)); - fprintf('The %.1f%% confidence interval for mean of h2>Q3 is (%.2f, %.2f). \n \n', 100*(1-alpha), CIs_h2(6,:)); - - fprintf('The %.1f%% confidence interval for mean of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(1,:)); - fprintf('The %.1f%% confidence interval for weighted mean of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(2,:)); - fprintf('The %.1f%% confidence interval for median (Q2) of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(3,:)); - fprintf('The %.1f%% confidence interval for the third quantile (Q3) of e2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(4,:)); - fprintf('The %.1f%% confidence interval for mean of e2>Q2 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(5,:)); - fprintf('The %.1f%% confidence interval for mean of e2>Q3 is (%.2f, %.2f). \n', 100*(1-alpha), CIs_e2(6,:)); - - save(fullfile(ACEfit_Par.ResDir,'Boot_CIs'),'alpha','CIs_h2','CIs_e2'); + fprintf('The %.1f%% confidence interval for mean of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(1,:)); + fprintf('The %.1f%% confidence interval for weighted mean of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(2,:)); + fprintf('The %.1f%% confidence interval for median (Q2) of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(3,:)); + fprintf('The %.1f%% confidence interval for the third quantile (Q3) of h2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(4,:)); + fprintf('The %.1f%% confidence interval for mean of h2>Q2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_h2(5,:)); + fprintf('The %.1f%% confidence interval for mean of h2>Q3 is (%.2f, %.2f). \n \n', 100*(1-Balpha), CIs_h2(6,:)); + + fprintf('The %.1f%% confidence interval for mean of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(1,:)); + fprintf('The %.1f%% confidence interval for weighted mean of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(2,:)); + fprintf('The %.1f%% confidence interval for median (Q2) of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(3,:)); + fprintf('The %.1f%% confidence interval for the third quantile (Q3) of e2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(4,:)); + fprintf('The %.1f%% confidence interval for mean of e2>Q2 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(5,:)); + fprintf('The %.1f%% confidence interval for mean of e2>Q3 is (%.2f, %.2f). \n', 100*(1-Balpha), CIs_e2(6,:)); + + save(fullfile(ACEfit_Par.ResDir,'Boot_CIs'),'Balpha','CIs_h2','CIs_e2'); end diff --git a/apace/CheckData.m b/apace/CheckData.m index 6bbcce1..eb70955 100644 --- a/apace/CheckData.m +++ b/apace/CheckData.m @@ -51,9 +51,10 @@ end % Number of all subjects -n = size(ACEfit_Par.kin,1); +n = ACEfit_Par.N; Vs.n = n; % redundant, but useful + % % Check family data and load the data as well (because you have to load % it to check it). diff --git a/apace/PrepParallel.m b/apace/PrepParallel.m index 5042f5c..05f407f 100644 --- a/apace/PrepParallel.m +++ b/apace/PrepParallel.m @@ -1,19 +1,23 @@ -function PrepParallel(ACEfit_Par,nParallel) +function PrepParallel(ACEfit_Par) % % Update the input structure array 'ACEfit_Par' for parallelization % -% nParallel - Number of parallel runs -% -if isempty(nParallel) - nParallel = 1; +if ( ~isfield(ACEfit_Par,'nPerm') || isempty(ACEfit_Par.nPerm) ) + ACEfit_Par.nPerm = 1000; +end +if ( ~isfield(ACEfit_Par,'nBoot') || isempty(ACEfit_Par.nBoot) ) + ACEfit_Par.nBoot = 1000; end +if ( ~isfield(ACEfit_Par,'nParallel') || isempty(ACEfit_Par.nParallel) ) + ACEfit_Par.nParallel = 1; +end +nPerm = ACEfit_Par.nPerm; +nBoot = ACEfit_Par.nBoot; +nParallel = ACEfit_Par.nParallel; ACEfit_Par.RunID = 1:nParallel; -nPerm = ACEfit_Par.nPerm; -nBoot = ACEfit_Par.nBoot; - nMZF = ACEfit_Par.nMZF; nDZF = ACEfit_Par.nDZF; nSibF = ACEfit_Par.nSibF; diff --git a/apace/ProcData.m b/apace/ProcData.m index 3cdce92..ba19ad1 100644 --- a/apace/ProcData.m +++ b/apace/ProcData.m @@ -3,14 +3,13 @@ % Process data for further computation % -[nElm,n] = size(Y); - -if ~isfield(ACEfit_Par,'Subset') || isempty(ACEfit_Par.Subset) - ACEfit_Par.Subset = 1:n; -end +nElm = size(Y,1); +n = ACEfit_Par.n; Subset = ACEfit_Par.Subset; +Y = Y(:,Subset); + X = ACEfit_Par.Dsnmtx; if isempty(X) X = ones(n,1); @@ -73,6 +72,11 @@ end end I_data = find(tmp); +if isempty(YM) && ~all(tmp==1) + fprintf('Voxels/elements removed due to missingness or being constant (%d).\n',sum(~tmp)); +elseif ~isempty(YM) && sum(~tmp)>sum(~YM) + fprintf('Voxels/elements removed due to missingness or being constant (%d; plus %d due to mask).\n',sum(~tmp)-sum(~YM),sum(~YM)); +end if isempty(I_data) error('All in-mask data vectors are null vectors or contain NaN') diff --git a/apace/Reorder.m b/apace/Reorder.m index 45e7a59..4e170ab 100644 --- a/apace/Reorder.m +++ b/apace/Reorder.m @@ -15,6 +15,11 @@ % 3 - Siblings % +% +% Fill-in defaults +% +if ~isfield(ACEfit_Par,'Subset'); ACEfit_Par.Subset = []; end + % % Read mother ID, father ID and zygosity % @@ -25,17 +30,55 @@ Afid = textscan(fid,'%s','Delimiter', ','); fclose(fid); -n = length(Afid{1})/4 - 1; -Bfid = reshape(Afid{1},4,n+1)'; +N = length(Afid{1})/4 - 1; +Bfid = reshape(Afid{1},4,N+1)'; Bfid = Bfid(2:end,:); -Bfid = cellfun(@(x) x(x~=' ' & x~='"'),Bfid,'un',0); +Bfid = cellfun(@(x) x(x~=' ' & x~='"'),Bfid,'UniformOutput',0); + +% Select a subset of subjects to form the sample for analysis +if ~isempty(ACEfit_Par.Subset) + Bfid = Bfid(ACEfit_Par.Subset,:); +else + ACEfit_Par.Subset = 1:N; +end + +% Save the total number of subjects in the kinship file +ACEfit_Par.N = N; + +% Compute and save the sample size after subject selection +n = length(ACEfit_Par.Subset); +ACEfit_Par.n = n; -FamIDs = cell2mat(cellfun(@str2num,strcat(Bfid(:,2),Bfid(:,3)),'un',0)); +%%% TODO: Error checking on Bfid to catch missing (i.e. empty or NaN) +%%% values for any input value, especially Father ID. +% Find missing data in the kinship file +Cfid = strcmp('NaN',Bfid) | strcmp('',Bfid); +if any(any(Cfid)) + fprintf('WARNING: Missing data (on %d subjects) found in the kinship file!\n',sum(any(Cfid,2)); +end + +% +% Construct Family ID using Mother ID and Father ID +% +FamIDs = cellfun(@str2num,strcat(Bfid(:,2),Bfid(:,3)),'UniformOutput',0); +FamIDs = cell2mat(cellfun(@(x) nansum(x),FamIDs,'UniformOutput',0)); +% Obtain subjects with missing Mother ID or Father ID +NaNSub = Cfid(:,2)+Cfid(:,3)~=0; +% Set Family ID of subjects with missing Mother ID or Father ID as +% max(FamIDs) + 1:length(NaNSub), i.e. subjects with missing kinship +% information are assumed to be singletons from different families +FamIDs(NaNSub) = max(FamIDs) + (1:sum(NaNSub))'; + +% +% Generate zygosity vector +% Zyg = zeros(n,1); Zyg(strcmp('MZ', Bfid(:,4))) = 1; Zyg(strcmp('NotMZ', Bfid(:,4))) = 2; Zyg(strcmp('NotTwin',Bfid(:,4))) = 3; +% Set subjects with missing Zygosity as singletons +Zyg(Zyg==0) = 3; % @@ -161,7 +204,7 @@ ' %4d Singletons\n',... ' %4d Smallest family size\n',... ' %4d Largest family size\n\n'],... - sum(nFam),nMZF,nDZF,nSibF,nSG,min(nFam),max(nFam)); + n,nMZF,nDZF,nSibF,nSG,min(nFam),max(nFam)); return diff --git a/apace/WK_ParallelOutline.m b/apace/WK_ParallelOutline.m index 78b2b17..6b118c4 100644 --- a/apace/WK_ParallelOutline.m +++ b/apace/WK_ParallelOutline.m @@ -8,104 +8,128 @@ % set(0, 'DefaultFigureVisible', 'off'); %%% 1) Create a structure array specifying details of data -ACEfit_Par.Model = 'ACE'; % Choose a model (AE or ACE) for - % data fitting. +ACEfit_Par.Model = 'ACE'; % Choose a model (AE or ACE) for + % data fitting. -ACEfit_Par.P_nm = 'APACE_imgs.txt'; % A list of images or other data - % specification; see FileFormats.txt +ACEfit_Par.P_nm = 'APACE_imgs.txt'; % A list of images or other data + % specification; see + % FileFormats.txt -ACEfit_Par.InfMx = 'KinInf.csv'; % Kinship information matrix of 4 - % columns with headers: - % 'SubjectID', 'MotherID', - % 'FatherID', and 'Zygosity'. - % MotherID and FatherID must be - % numeric, and Zygosity is 'MZ', - % 'NotMZ' or 'NotTwin'. The order - % of subjects must match that of - % image paths or subjects in - % ACEfit_Par.P_nm. -ACEfit_Par.ResDir = '/my/path/ResDir'; +ACEfit_Par.InfMx = 'KinInf.csv'; % Kinship information CSV file + % with header row and exactly 4 + % columns: Subject ID, Mother ID, + % Father ID, and Zygosity. Mother + % ID and Father ID must be + % positive integers, and Zygosity + % is 'MZ', 'NotMZ' or 'NotTwin'. + % The order of subjects must + % match that of image paths or + % subjects in ACEfit_Par.P_nm. + +ACEfit_Par.ResDir = './ResDir'; % Set result folder %%% The rest are optional; omit to use default values -ACEfit_Par.Subset = []; % Set equal to (subset of) - % subject indicies to be - % considered in the - % analysis. Empty (default) means - % "all subjects". - -ACEfit_Par.Pmask = 'APACE_mask.nii'; % Brain mask image (default: whole volume) - -ACEfit_Par.Dsnmtx = ''; % Design matrix (default: all-ones vector) - -ACEfit_Par.Nlz = 1; % Inverse Gaussian normalisation options - % (applied to each phenotype voxel/element) - % 0 - None - % 1 - Gaussian normalisation *before* - % forming residuals - % 2 - Gaussian normalisation *after* - % forming residuals; note, though, - % that after this Gaussian - % normalisation the residuals - % will be formed again to ensure - % orthognality of residuals to - % the design matrix. - -ACEfit_Par.AggNlz = 0; % Aggregate heritability normalisation options - % (applied to each phenotype voxel/element) - % 0 - Default, meaning that only - % determined by Nlz option. - % This always entails - % de-meaning and removal of any effects - % in Dsnmtx, but variance at each - % voxel/element unchanged. - % 1 - Same as 0, but with variance - % normalisation; this is full - % 'studentization'. - % 2 - *Undo* mean centering of - % residual formation. (If Dsnmtx - % is empty or default value of - % an all-ones vector, this is - % equivalent to using the raw - % input data. If Dsnmtx - % contains nuisance regressors - % the residuals will be formed - % and the mean added back in. - % 3 - Same as 2, but with - % variance normalisation; - % note that for each - % voxel/element, mean is - % first added back and then - % data are divided by stdev. - -ACEfit_Par.ContSel = []; % Select a single contrast (a volume - % in a 4D Nifti file supplied for - % each subject, or at the last - % dimension in a cifti image; NOT - % compatible with a single file - % containing all subjects' data.) +ACEfit_Par.Subset = []; % Set equal to (a subset of) + % subject indicies to be + % considered in the analysis. + % Empty (default) means "all + % subjects". + +ACEfit_Par.Pmask = ''; % Brain mask image (default: + % whole volume) + +ACEfit_Par.Dsnmtx = ''; % Design matrix (default: + % all-ones vector) + +ACEfit_Par.Nlz = 1; % Inverse Gaussian normalisation + % options (applied to each + % phenotype voxel/element) + % 0 - None + % 1 - Gaussian normalisation + % *before* forming residuals + % (default) + % 2 - Gaussian normalisation + % *after* forming residuals; + % note, though, that after + % this Gaussian + % normalisation the + % residuals will be formed + % again to ensure + % orthognality of residuals + % to the design matrix. + +ACEfit_Par.AggNlz = 0; % Aggregate heritability + % normalisation options (applied + % to each phenotype + % voxel/element) + % 0 - Default, meaning that only + % determined by Nlz option. + % This always entails + % de-meaning and removal of + % any effects in Dsnmtx, but + % variance at each + % voxel/element unchanged. + % 1 - Same as 0, but with + % variance normalisation; + % this is full + % 'studentization'. + % 2 - *Undo* mean centering of + % residual formation. (If + % Dsnmtx is empty or default + % value of an all-ones + % vector, this is equivalent + % to using the raw input + % data. If Dsnmtx contains + % nuisance regressors the + % residuals will be formed + % and the mean added back + % in. + % 3 - Same as 2, but with + % variance normalisation; + % note that for each + % voxel/element, mean is + % first added back and then + % data are divided by stdev. + +ACEfit_Par.ContSel = []; % Select a single contrast (a + % volume in a 4D Nifti file + % supplied for each subject, or + % at the last dimension in a + % cifti image; NOT compatible + % with a single file containing + % all subjects' data.) -ACEfit_Par.NoImg = 0; % If 1, suppress image-wise inference, - % and only compute summaries. +ACEfit_Par.NoImg = 0; % If 1, suppress image-wise + % inference, and only compute + % summaries. +ACEfit_Par.alpha_CFT = []; % Cluster-forming threshold + % (default: 0.05) + +ACEfit_Par.nPerm = []; % Number of permutations + % (default: 1000) + +ACEfit_Par.nBoot = []; % Number of bootstrap replicates + % (default: 1000) + +ACEfit_Par.nParallel = []; % Number of parallel runs + % (default: 1, without + % parallelization) -%%% 2) Updata 'ACEfit_Par' with the input data information +%%% 2) Update 'ACEfit_Par' with the input data information ACEfit_Par = PrepData(ACEfit_Par); %%% 3) Run the original data once -ACEfit_Par.alpha_CFT = []; % Cluster-forming threshold (default: 0.05) ACEfit_Par = ACEfit(ACEfit_Par); -%%% 4) Add permutation and bootstrapping information, and save "ACEfit_Par.mat" -ACEfit_Par.nPerm = 1000; % Number of permutations -ACEfit_Par.nBoot = 1000; % Number of bootstrap replicates -nParallel = []; % Number of parallel runs (default: 1, without parallelization) - -PrepParallel(ACEfit_Par,nParallel); +%%% 4) Add permutation and bootstrapping information, and save +%%% "ACEfit_Par.mat" +PrepParallel(ACEfit_Par); % -% Permutation inference for computing FWE-corrected p-values +% Permutation inference for computing FWE/FDR-corrected p-values % if ACEfit_Par.nPerm>0 @@ -122,15 +146,14 @@ % one for each RunID. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%% 2) This will merge together all the results from the specified RunID's + %%% 2) This will merge together all the results from the specified + %%% RunID's load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat')); ACEfit_Perm_Parallel_Results(ACEfit_Par); %%% 3) - FWEalpha = 0.05; - FDRalpha = 0.05; load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat')); - ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha); + ACEfit_Results(ACEfit_Par); end @@ -157,10 +180,9 @@ load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat')); ACEfit_Boot_Parallel_Results(ACEfit_Par); - %%% 3) - Balpha = 0.05; + %%% 3) Construct CIs load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat')); - Boot_CIs(ACEfit_Par,Balpha) + Boot_CIs(ACEfit_Par); end @@ -173,24 +195,24 @@ % skipped by setting ACEfit_Par.nPerm=0 and ACEfit_Par.nBoot=0 separately; % the following code can be run immediately after "PrepParallel". % + load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat')); -Palpha = 0.05; % Significance threhold for permutations (plotting only) -Balpha = 0.05; % Confidence level, where CI's have level 100*(1-alpha) +AgHe_Method(ACEfit_Par); -AgHe_Method(ACEfit_Par,Palpha,Balpha) % % Once with no variance normalisation -% ACEfit_Par.AggNlz = 0; % de-meaning only -% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_NoNorm') +% ACEfit_Par.AggNlz = 0; % de-meaning only +% AgHe_Method(ACEfit_Par,'_NoNorm'); % -% % Now, again, but with Variance normalisation -% ACEfit_Par.AggNlz = 1; % de-meaning and scaling to have stdev of 1.0 -% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_Norm') +% % Now, again, but with variance normalisation +% ACEfit_Par.AggNlz = 1; % de-meaning and scaling to have stdev of 1.0 +% AgHe_Method(ACEfit_Par,'_Norm'); % % Generate summary file % -APACEsummary(ACEfit_Par,'ResultSummary') +APACEsummary(ACEfit_Par,'ResultSummary'); % Save results to "ResultsSummary.csv" +% APACEsummary(ACEfit_Par); % Print results to the screen