lab:v1_on_nov_23_2008
/*********************************************************** *Mediation Analysis * * * * Author: Zhiyong Zhang * * * * Purpose: Test mediation through bootstrap and EM * * algorithm for missing data * * * * Usage: * * Run the program and the message appears * ***********************************************************/ /*********************************************************** Change history: 2008/12/20 Version 1 based on the BootMed **********************************************************/ // The library newmat was used for the matrix operation // #define WANT_STREAM #define WANT_MATH #define WANT_FSTREAM #include "newmatap.h" #include "newmatio.h" #include "newmatnl.h" #include "newran.h" #include "math.h" #include<iostream> #include<fstream> #include<ostream> #include<string> #include "time.h" #ifdef use_namespace using namespace RBD_LIBRARIES; #endif /*------------------------------------------------------*\ * Normal quantile calculation function * \*------------------------------------------------------*/ double dfQNorm(double alpha){ ColumnVector ck(100); ck=0.0; ck(1)=1; ck(2)=1; for (int k=3; k<=100; k++){ for (int m=1; m<k; m++){ ck(k)=ck(k)+ck(m)*ck(k-m)/((m)*(2.0*m-1.0)); } } double x; x=2.0*alpha-1.0; double pi=3.14159265358979323846; double erf=0.0; for (int i=1; i<=100; i++){ erf=erf+ck(i)/(2.0*i-1.0)*pow((sqrt(pi)/2.0*x),(2.0*i-1.0)); } erf=sqrt(2.0)*erf; return erf; } // factorial fucntion double ifFactor(int N){ double iFactor=1.0; for (int i=1; i<=N; i++){ iFactor = iFactor * i; } return iFactor; } /*------------------------------------------------------*\ * Normal cumulative distribution function * \*------------------------------------------------------*/ double dfPNorm(double dX){ double dZ = dX/sqrt(2.0); double pi=3.14159265358979323846; double erf=dZ; for (int i=1; i<=100; i++){ erf=erf+pow(-1.0,i)*pow(dZ, 2.0*i+1.0)/(ifFactor(i)*(2.0*i+1.0)); } erf=2.0*erf/sqrt(pi); erf=(1.0+erf)/2.0; return erf; } /*------------------------------------------------------*\ * Function to read in a data file and decide the size. * * Read in the file name and output the row and column * * Numbers. * \*------------------------------------------------------*/ void size(char *filename, int &rn, int &cn) { Tracer tr("size"); // for tracking exceptions //Open file and decide the size of the file std::ifstream in; in.open(filename); if (!in) { std::cout << "cannot open the data file\n"; exit(1); } std::string str; std::string buf; // Decide the number of rows int row=0; while(std::getline(in,str)) { if(row==0) buf=str; row++; // std::cout<<str<<"\n"; } // Number of columns //c_str() convert a string to char* int column=0; const char *c=buf.c_str(); bool tag=true; while(*c!='\0') { if(*c!=' '&&tag==true) { c++; tag=false; column++; } else if(*c++==' ') { tag=true; } } rn=row; cn=column; in.close(); } /*------------------------------------------------------*\ * Read the data from a file into a matrix * \*------------------------------------------------------*/ Matrix readdata(char *filename, int row, int col){ std::ifstream fin; fin.open(filename); if (!fin) { std::cout << "cannot open the data file\n"; exit(1); } Matrix X(row,col); for (int i=1; i<=row; i++) for (int j=1; j<=col; j++) { fin >> X(i,j); } fin.close(); return X; } /*------------------------------------------------------*\ * Function to organize the data using missing pattern * * Return the matrix of data and vector of pattern size * * Input: matData: original data * * iN: sample size * * Output: cvN: a vector with pattern size * * matPatData: sorted data with missing pattern * \*------------------------------------------------------*/ void pattern(Matrix matData, int iN, ColumnVector &cvN, Matrix &matPatData){ int n1=0,n2=0,n3=0,n4=0,n5=0,n6=0,n7=0,n8=0; double dmiss=99999.0; //Get the dimension of different patterns for (int i1=1; i1<=iN; i1++){ if (matData(i1,1)!=dmiss){ if (matData(i1,2)!=dmiss){ if (matData(i1,3)!=dmiss){ n1=n1+1; }else{ n2=n2+1; } }else{ if (matData(i1,3)!=dmiss){ n3=n3+1; }else{ n4=n4+1; } } }else{ if (matData(i1,2)!=dmiss){ if (matData(i1,3)!=dmiss){ n5=n5+1; }else{ n6=n6+1; } }else{ if (matData(i1,3)!=dmiss){ n7=n7+1; }else{ n8=n8+1; } } } }//end of the sample size for each pattern //Define the matrix to store the data for different patterns int m1=n1, m2=n2, m3=n3, m4=n4, m5=n5, m6=n6, m7=n7, m8=n8; if (n1==0){m1=1;} if (n2==0){m2=1;} if (n3==0){m3=1;} if (n4==0){m4=1;} if (n5==0){m5=1;} if (n6==0){m6=1;} if (n7==0){m7=1;} if (n8==0){m8=1;} Matrix matPat1(iN,3); Matrix matPat2(iN,3); Matrix matPat3(iN,3); Matrix matPat4(iN,3); Matrix matPat5(iN,3); Matrix matPat6(iN,3); Matrix matPat7(iN,3); Matrix matPat8(iN,3); matPat1=0.0; matPat2=0.0; matPat3=0.0; matPat4=0.0; matPat5=0.0; matPat6=0.0; matPat7=0.0; matPat8=0.0; //Get the dimension of different patterns int k1=1, k2=1, k3=1, k4=1, k5=1, k6=1, k7=1, k8=1; for (int i2=1; i2<=iN; i2++){ if (matData(i2,1)!=dmiss){ if (matData(i2,2)!=dmiss){ if (matData(i2,3)!=dmiss){ matPat1.row(k1)=matData.row(i2); k1++; }else{ matPat2.row(k2)=matData.row(i2); k2++; } }else{ if (matData(i2,3)!=dmiss){ matPat3.row(k3)=matData.row(i2); k3++; }else{ matPat4.row(k4)=matData.row(i2); k4++; } } }else{ if (matData(i2,2)!=dmiss){ if (matData(i2,3)!=dmiss){ matPat5.row(k5)=matData.row(i2); k5++; }else{ matPat6.row(k6)=matData.row(i2); k6++; } }else{ if (matData(i2,3)!=dmiss){ matPat7.row(k7)=matData.row(i2); k7++; }else{ matPat8.row(k8)=matData.row(i2); k8++; } } } }//end of the data pattern //put the data into the pattern size vector cvN<<n1<<n2<<n3<<n4<<n5<<n6<<n7<<n8; //put the data into the pattern matrix int l2=n1+n2; int l3=l2+n3; int l4=l3+n4; int l5=l4+n5; int l6=l5+n6; int l7=l6+n7; int i3; if (n1>0){ for (i3=1; i3<=n1; i3++){ matPatData.row(i3)=matPat1.row(i3); } } if (n2>0){ for (i3=1; i3<=n2; i3++){ matPatData.row(n1+i3)=matPat2.row(i3); } } if (n3>0){ for (i3=1; i3<=n3; i3++){ matPatData.row(l2+i3)=matPat3.row(i3); } } if (n4>0){ for (i3=1; i3<=n4; i3++){ matPatData.row(l3+i3)=matPat4.row(i3); } } if (n5>0){ for (i3=1; i3<=n5; i3++){ matPatData.row(l4+i3)=matPat5.row(i3); } } if (n6>0){ for (i3=1; i3<=n6; i3++){ matPatData.row(l5+i3)=matPat6.row(i3); } } if (n7>0){ for (i3=1; i3<=n7; i3++){ matPatData.row(l6+i3)=matPat7.row(i3); } } if (n8>0){ for (i3=1; i3<=n8; i3++){ matPatData.row(l7+i3)=matPat8.row(i3); } } } /*------------------------------------------------------*\ * Function to calculate the mean and covariance matrix * * Output: mCov mCov[1:c,r] is the the means mCov[1:c,1:r-1] is the covariance \*------------------------------------------------------*/ Matrix mfCov(Matrix matX, int iRow, int iCol){ RowVector vMean = matX.sum_columns()/iRow; Matrix mXbar(iRow,iCol); for (int i=1; i<=iRow; i++){ mXbar.row(i)=vMean; } Matrix mDesign=matX-mXbar; Matrix mCov(iCol+1, iCol); mCov.submatrix(1,iCol,1,iCol)=mDesign.t()*mDesign/(iRow-1.0); mCov.row(iCol+1)=vMean; return mCov; } /*------------------------------------------------------*\ * Function to estimate the coefficients * \*------------------------------------------------------*/ RowVector cvfParas(Matrix mCov) { double d_a1, d_a2, d_b1, d_sig1, d_sig2, d_sigx, d_a2b1, dybar, dmbar, dxbar; double d_sy, d_sx, d_sm, d_sxy, d_sxm, d_sym; double d_mM, d_mY; RowVector cvPara(8); d_sy=mCov(3,3); d_sx=mCov(1,1); d_sm=mCov(2,2); d_sxy=mCov(1,3); d_sym=mCov(2,3); d_sxm=mCov(1,2); dybar=mCov(4,3); dxbar=mCov(4,1); dmbar=mCov(4,2); if (d_sxm*d_sxm-d_sm*d_sx==0){ cout<<"Problem in calculation"<<endl; exit(1); }else{ d_a1=(d_sxm*d_sym-d_sm*d_sxy)/(d_sxm*d_sxm-d_sm*d_sx); d_a2=(d_sx*d_sym-d_sxm*d_sxy)/(d_sm*d_sx-d_sxm*d_sxm); d_b1=d_sxm/d_sx; d_a2b1=d_a2*d_b1; d_sigx=d_sx; d_sig1=d_sm-d_sxm*d_sxm/d_sx; d_sig2=(d_sm*d_sx*d_sy-d_sxm*d_sxm*d_sy-d_sx*d_sym*d_sym-d_sm*d_sxy*d_sxy+2*d_sxm*d_sym*d_sxy)/(d_sm*d_sx-d_sxm*d_sxm); d_mM = dmbar - dxbar*d_b1; d_mY = dybar - d_a1*dxbar - d_a2*dmbar; cvPara(1)=d_mM; cvPara(2)=d_mY; cvPara(3)=d_b1; cvPara(4)=d_a2; cvPara(5)=d_a1; //variance for residuals of M~X cvPara(6)=d_sig1; //variance for residuals of Y~X+M cvPara(7)=d_sig2; cvPara(8)=d_a2b1; } return cvPara; } /*------------------------------------------------------*\ * EM algorithm - The M-step * Output: mCov -- covariance matrix and means \*------------------------------------------------------*/ Matrix mfMstep(Matrix mSS, int iN){ Matrix mSSV(3,3); RowVector mSSM(3); mSSV = mSS.submatrix(1,3,1,3); mSSM = mSS.row(4); RowVector mMean=mSSM/iN; Matrix mCov(4,3); mCov.submatrix(1,3,1,3) = mSSV/(iN-1.0) - iN/(iN-1.0)*mMean.t()*mMean; mCov.row(4)=mMean; return mCov; } /*------------------------------------------------------*\ * EM algorithm - The E-step * Output: mCov -- covariance matrix and means \*------------------------------------------------------*/ Matrix mfEstep(Matrix matData, Matrix mCov, int iN){ ColumnVector sX2(iN), sM2(iN),sY2(iN),sXM(iN),sXY(iN),sMY(iN),mX(iN),mM(iN),mY(iN); int i; double dmiss=99999.0; //Get the dimension of different patterns for (i=1; i<=iN; i++){ if (matData(i,1)!=dmiss){ if (matData(i,2)!=dmiss){ if (matData(i,3)!=dmiss){ mX(i)=matData(i,1); mM(i)=matData(i,2); mY(i)=matData(i,3); sX2(i)=mX(i)*mX(i); sY2(i)=mY(i)*mY(i); sM2(i)=mM(i)*mM(i); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i); }else{ mX(i)=matData(i,1); mM(i)=matData(i,2); mY(i)=mCov(4,3) + ( (mCov(3,1)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mX(i)-mCov(4,1))+ (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*(mM(i)-mCov(4,2)))/( mCov(1,1)*mCov(2,2)-mCov(1,2)*mCov(1,2)); sX2(i)=mX(i)*mX(i); sY2(i)=mY(i)*mY(i)+ mCov(3,3)-( (mCov(3,1)*mCov(2,2)-mCov(3,2)*mCov(1,2))*mCov(3,1)+ (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*mCov(3,2))/( mCov(1,1)*mCov(2,2)-mCov(1,2)*mCov(1,2)); sM2(i)=mM(i)*mM(i); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i); } }else{ if (matData(i,3)!=dmiss){ mX(i)=matData(i,1); mY(i)=matData(i,3); mM(i)=mCov(4,2)+( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mX(i)-mCov(4,1))+ (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*(mY(i)-mCov(4,3)))/( mCov(1,1)*mCov(3,3)-mCov(1,3)*mCov(1,3)); sX2(i)=mX(i)*mX(i); sY2(i)=mY(i)*mY(i); sM2(i)=mM(i)*mM(i)+mCov(2,2)-( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*mCov(2,1)+ (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*mCov(2,3))/( mCov(1,1)*mCov(3,3)-mCov(1,3)*mCov(1,3)); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i); }else{ mX(i)=matData(i,1); mM(i)=mCov(4,2)+mCov(1,2)/mCov(1,1)*(mX(i)-mCov(4,1)); mY(i)=mCov(4,3)+mCov(1,3)/mCov(1,1)*(mX(i)-mCov(4,1)); sX2(i)=mX(i)*mX(i); sY2(i)=mY(i)*mY(i)+mCov(3,3)-mCov(1,3)*mCov(1,3)/mCov(1,1); sM2(i)=mM(i)*mM(i)+mCov(2,2)-mCov(1,2)*mCov(1,2)/mCov(1,1); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i)+mCov(2,3)-mCov(1,2)*mCov(1,3)/mCov(1,1); } } }else{ if (matData(i,2)!=dmiss){ if (matData(i,3)!=dmiss){ mM(i)=matData(i,2); mY(i)=matData(i,3); mX(i)=mCov(4,1) + ( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mM(i)-mCov(4,2))+ (mCov(1,3)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mY(i)-mCov(4,3)))/( mCov(2,2)*mCov(3,3)-mCov(2,3)*mCov(2,3)); sX2(i)=mX(i)*mX(i) + mCov(1,1)-( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mCov(1,2))+ (mCov(1,3)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mCov(1,3)))/( mCov(2,2)*mCov(3,3)-mCov(2,3)*mCov(2,3)); sY2(i)=mY(i)*mY(i); sM2(i)=mM(i)*mM(i); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i); }else{ mM(i)=matData(i,2); mX(i)=mCov(4,1)+mCov(1,2)/mCov(2,2)*(mM(i)-mCov(4,2)); mY(i)=mCov(4,3)+mCov(2,3)/mCov(2,2)*(mM(i)-mCov(4,2)); sX2(i)=mX(i)*mX(i)+mCov(1,1)-mCov(1,2)*mCov(1,2)/mCov(2,2); sY2(i)=mY(i)*mY(i)+mCov(3,3)-mCov(3,2)*mCov(3,2)/mCov(2,2); sM2(i)=mM(i)*mM(i); sXM(i)=mX(i)*mM(i); sXY(i)=mX(i)*mY(i)+mCov(1,3)-mCov(1,2)*mCov(2,3)/mCov(2,2); sMY(i)=mM(i)*mY(i); } }else{ if (matData(i,3)!=dmiss){ mY(i)=matData(i,3); mX(i)=mCov(4,1)+mCov(1,3)/mCov(3,3)*(mY(i)-mCov(4,3)); mM(i)=mCov(4,2)+mCov(2,3)/mCov(3,3)*(mY(i)-mCov(4,3)); sX2(i)=mX(i)*mX(i)+mCov(1,1)-mCov(1,3)*mCov(1,3)/mCov(3,3); sY2(i)=mY(i)*mY(i); sM2(i)=mM(i)*mM(i)+mCov(2,2)-mCov(2,3)*mCov(2,3)/mCov(3,3); sXM(i)=mX(i)*mM(i)+mCov(1,2)-mCov(1,3)*mCov(2,3)/mCov(3,3); sXY(i)=mX(i)*mY(i); sMY(i)=mM(i)*mY(i); } } } }//end of the sample size for each pattern //Get the sum of each vector and put them into a matrix Matrix mSS(4,3); mSS(1,1)=sX2.sum(); mSS(1,2)=sXM.sum(); mSS(1,3)=sXY.sum(); mSS(2,1)=mSS(1,2); mSS(2,2)=sM2.sum(); mSS(2,3)=sMY.sum(); mSS(3,1)=mSS(1,3); mSS(3,2)=mSS(2,3); mSS(3,3)=sY2.sum(); mSS(4,1)=mX.sum(); mSS(4,2)=mM.sum(); mSS(4,3)=mY.sum(); return mSS; } // Function to calcualte the difference of two matrix double dfDiff(Matrix matCov1, Matrix matCov2){ Matrix mDiff(4,3); double dError; int i,j; for (i=1; i<=4; i++){ for (j=1; j<=3; j++){ mDiff(i,j)=(matCov1(i,j)-matCov2(i,j))/matCov2(i,j); } } dError=mDiff.sum_absolute_value(); return dError; } /*------------------------------------------------------*\ * EM algorithm * \*------------------------------------------------------*/ RowVector cvfEM(Matrix matData, Matrix mCov, int iN){ RowVector cvPara; double e=1.0; Matrix ssCov(4,3); Matrix mTCov(4,3); mTCov=mCov; Matrix mT2Cov(4,3); int iMaxI=10000; int j=1; while (e>0.0001 || j==iMaxI){ ssCov=mfEstep(matData, mTCov, iN); mT2Cov=mfMstep(ssCov,iN); e=dfDiff(mT2Cov, mTCov); mTCov=mT2Cov; j++; } if (j==iMaxI){ cout<<"\nThe EM algorithm cannot converge..."<<endl; exit(1); } cvPara=cvfParas(mTCov); return cvPara; } /*------------------------------------------------------*\ * Delete a row from a matrix * \*------------------------------------------------------*/ Matrix mDelARow(Matrix matData, int iR, int iN){ Matrix mDel(iN-1, 3); if (iR==1){ mDel=matData.submatrix(2,iN,1,3); }else if (iR==iN){ mDel=matData.submatrix(1,iN-1,1,3); }else{ mDel.submatrix(1,iR-1,1,3)=matData.submatrix(1,iR-1,1,3); mDel.submatrix(iR, iN-1,1,3)=matData.submatrix(iR+1,iN,1,3); } return mDel; } /*------------------------------------------------------*\ * Calculate the acceleration factor using Jackknife * \*------------------------------------------------------*/ RowVector rvJackA(Matrix matData, Matrix mCov,int iN){ Matrix mJackPara(iN, 8), mDel(iN-1,3), mTemp1(iN,8), mTemp2(iN,8); RowVector dA(8), dMean(8), dSum1(8), dSum2(8); int i,j; for (i=1;i<=iN;i++){ mDel = mDelARow(matData,i,iN); mJackPara.row(i)=cvfEM(mDel,mCov,iN-1); } dMean = mJackPara.sum_columns()/iN; for (i=1; i<=iN; i++){ for (j=1; j<=8;j++){ mTemp1(i,j) = (dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j)); mTemp2(i,j) = (dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j)); } } dSum1=mTemp1.sum_columns(); dSum2=mTemp2.sum_columns(); for (i=1; i<=8; i++){ dA(i)=dSum1(i)/(pow(dSum2(i),1.5)*6.0); } return dA; } /*------------------------------------------------------*\ * Function to bootstrap * \*------------------------------------------------------*/ Matrix mfBootOnce(Matrix matData, int iRow){ Matrix mBoot(iRow, 3); Uniform U; int intvalue; for (int i=1; i<=iRow; i++){ intvalue=int(U.Next()*iRow)+1; mBoot.row(i)=matData.row(intvalue); } return mBoot; } /*------------------------------------------------------*\ * The stratified bootstrap * \*------------------------------------------------------*/ Matrix mfBootstrap(Matrix matData, ColumnVector cvN, int iN){ Matrix mStratBoot(iN,3); int n1=(int)cvN(1); int n2=(int)cvN(2); int n3=(int)cvN(3); int n4=(int)cvN(4); int n5=(int)cvN(5); int n6=(int)cvN(6); int n7=(int)cvN(7); int m2=n1+n2; int m3=m2+n3; int m4=m3+n4; int m5=m4+n5; int m6=m5+n6; int m7=m6+n7; if (n1>0){ mStratBoot.submatrix(1,n1,1,3)=mfBootOnce(matData.submatrix(1,n1,1,3),n1); } if (n2>0){ mStratBoot.submatrix(n1+1,m2,1,3)=mfBootOnce(matData.submatrix(n1+1,m2,1,3),n2); } if (n3>0){ mStratBoot.submatrix(m2+1,m3,1,3)=mfBootOnce(matData.submatrix(m2+1,m3,1,3),n3); } if (n4>0){ mStratBoot.submatrix(m3+1,m4,1,3)=mfBootOnce(matData.submatrix(m3+1,m4,1,3),n4); } if (n5>0){ mStratBoot.submatrix(m4+1,m5,1,3)=mfBootOnce(matData.submatrix(m4+1,m5,1,3),n5); } if (n6>0){ mStratBoot.submatrix(m5+1,m6,1,3)=mfBootOnce(matData.submatrix(m5+1,m6,1,3),n6); } if (n7>0){ mStratBoot.submatrix(m6+1,m7,1,3)=mfBootOnce(matData.submatrix(m6+1,m7,1,3),n7); } return mStratBoot; } /*------------------------------------------------------*\ * Function to calculate the CI for one parameter * \*------------------------------------------------------*/ RowVector rvfCI(ColumnVector rvBootRes, double dThetaHat, double dA, int iB, double alpha){ RowVector rvCI(6); int i, iUp, iLow; double iPos=0.0; double tt1, tt2; if (alpha < 1.0-alpha){ alpha = 1.0-alpha; } sort_ascending(rvBootRes); for (i=1; i<=iB; i++){ if (rvBootRes(i)<dThetaHat){ iPos=iPos+1.0; } } iLow = int((1.0-alpha)/2.0*iB+.5); rvCI(1)=rvBootRes(iLow); iUp = int((1.0+alpha)/2.0*iB+.5); rvCI(2)=rvBootRes(iUp); double dZ01=dfQNorm(iPos/iB); double zalpha1=dfQNorm((1.0-alpha)/2.0); double zalpha2=dfQNorm((1.0+alpha)/2.0); //CI for BC tt1=dfPNorm(2.0*dZ01+zalpha1); iLow = int(tt1*iB+.5); rvCI(3)=rvBootRes(iLow); tt2=dfPNorm(2.0*dZ01+zalpha2); iUp = int(tt2*iB+.5); rvCI(4)=rvBootRes(iUp); //CI for BCa tt1 = dfPNorm(dZ01+(dZ01+zalpha1)/(1-dA*(dZ01+zalpha1))); iLow = int(tt1*iB+.5); rvCI(5)=rvBootRes(iLow); tt2 = dfPNorm(dZ01+(dZ01+zalpha2)/(1-dA*(dZ01+zalpha2))); iUp = int(tt2*iB+.5); rvCI(6)=rvBootRes(iUp); return rvCI; } /*------------------------------------------------------*\ * Function to calculate the CIs * \*------------------------------------------------------*/ Matrix mBootCI(Matrix matData, Matrix mCov, ColumnVector cvN, int iN, int iB,double alpha){ RowVector rvTrue=cvfEM(matData,mCov,iN); RowVector rvTemp(6), tempres(8); Matrix mBootData(iN,3),mBootRes(iB,8), mCIres(8,7); int k=int(iB/10); int i,j=1; cout<<"\nRunning: "; for (i=1; i<=iB; i++){ mBootData=mfBootstrap(matData,cvN,iN); tempres=cvfEM(mBootData,mCov,iN); mBootRes.row(i)=tempres; if (i==j*k){ cout<<j*10<<"% "; j++; } } //Run Jackknife cout<<"Finishing "; RowVector rvA=rvJackA(matData,mCov,iN); cout<<"DONE"<<"\n"; for (i=1; i<=8; i++){ rvTemp = rvfCI(mBootRes.column(i),rvTrue(i),rvA(i),iB,alpha); mCIres(i,1)=rvTrue(i); mCIres(i,2)=rvTemp(1); mCIres(i,3)=rvTemp(2); mCIres(i,4)=rvTemp(3); mCIres(i,5)=rvTemp(4); mCIres(i,6)=rvTemp(5); mCIres(i,7)=rvTemp(6); } return mCIres; } /*------------------------------------------------------*\ * The main program * \*------------------------------------------------------*/ int main() { // Print the program information cout <<"\n"; cout <<"\n"; cout <<"---------------------------------------------------------"<<endl; cout <<"| Program name: BootMedEM.exe (V1.0) |"<<endl; cout <<"| By Zhiyong Zhang & Lijuan Wang |"<<endl; cout <<"| Email zzhang4@nd.edu for questions and comments |"<<endl; cout <<"| See readme.txt for more information |"<<endl; cout <<"---------------------------------------------------------"<<endl; cout <<"\n"; cout <<"\n"; cout <<"*********************************************************"<<endl; cout <<"* WARNING! *"<<endl; cout <<"* Although this program has been tested to perform as *"<<endl; cout <<"* expected in many cases, we cannot guarantee its *"<<endl; cout <<"* performance under all possible circumstances. *"<<endl; cout <<"* YOU CAN USE IT FOR FREE BUT AT YOUR OWN RISK. *"<<endl; cout <<"*********************************************************"<<endl; cout <<endl<<endl; // Ask the input of the output file cout <<"Please input the output file name: "<<endl; std::string output; std::cin >>output; char *fout; fout=(char*)output.c_str(); std::ofstream out; out.open(fout); // Print out the file information out <<"----------------------------------------------------------\n"; out <<"| Program name: BootMedEM.exe (V1.0) |\n"; out <<"| By Zhiyong Zhang & Lijuan Wang |\n"; out <<"| Email zzhang4@nd.edu for questions and comments |\n"; out <<"| See readme.txt for more information |\n"; out <<"----------------------------------------------------------\n"; out <<"\n"; out <<"\n"; out <<"**********************************************************\n"; out <<"* WARNING! *\n"; out <<"* Although this program has been tested to perform as *\n"; out <<"* expected in many cases, we cannot guarantee its *\n"; out <<"* performance under all possible circumstances. *\n"; out <<"* YOU CAN USE IT FOR FREE BUT AT YOUR OWN RISK. *\n"; out <<"**********************************************************\n"; out <<"\n\n"; //time time_t rawtime; time(&rawtime); out <<"The BootMed program is ran on "<<ctime(&rawtime)<<"\n"; out << "The output file is: "; out << output << "\n"; //Input the data file name and check the data int c1; int r1; std::string filename; cout <<"\nPlease input the file name of the data:"<<endl; out <<"The input the file of the data is: "; std::cin >> filename; out << filename <<"\n"; char *fstr; fstr= (char*)filename.c_str(); size(fstr,r1, c1); //Decide if number of variables > 3 if (c1>3){ cout<<"The number of variables >3 !"; cout<<"Only three varialbes (Y, X, Mediation) are allowed in the data file!"; out<<"The number of variables >3 !\n"; out<<"Only three varialbes (Y, X, Mediation) are allowed in the data file!"; exit(1); } //Read the data into a matrix Matrix INPORG=readdata(fstr,r1,c1); //request the input of the alpha level std::string stralpha; cout <<"\nPlease specify the alpha level (a number between 0 and 1 only):"<<endl; out <<"The alpha level is: "; std::cin >>stralpha; double alpha= atof(stralpha.c_str()); //cout <<alpha <<endl; out <<alpha<<"\n"; if (alpha>=1 || alpha<=0){ cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<<endl; out <<"Only a number bigger than 0 and less than 1 is accepted for the seed!\n"; exit(1); } //request the input of seed for random number generation std::string strseed; cout <<"\nPlease specify the random number seed (a number between 0 and 1 only):"<<endl; out <<"The random number seed is: "; std::cin >>strseed; out <<strseed<<"\n"; double seed= atof(strseed.c_str()); //cout <<seed; if (seed>=1 || seed<=0){ cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<<endl; out <<"Only a number bigger than 0 and less than 1 is accepted for the seed!\n"; exit(1); } MotherOfAll urng(seed); // declare uniform random number generator Random::Set(urng); //Calculate the patterns of missing data ColumnVector cvN(8); Matrix matPatData(r1,c1); pattern(INPORG,r1,cvN,matPatData); //output the missing data pattern int iP=1; out <<"\n"; out <<"--------------------------------"<<"\n"; out <<"| Missing Data Pattern |"<<"\n"; out <<"--------------------------------"<<"\n"; out <<"\n"; out.width(7);out <<"Pattern"; out.width(3);out<<"X"; out.width(3);out<<"M"; out.width(3);out<<"Y"<<"\n"; if (cvN(1)>0){ out.width(1);out<<iP;out.width(16);out<<" o o o\n"; iP++; } if (cvN(2)>0){ out.width(1);out<<iP;out.width(16);out<<" o o x\n"; iP++; } if (cvN(3)>0){ out.width(1);out<<iP;out.width(16);out<<" o x o\n"; iP++; } if (cvN(4)>0){ out.width(1);out<<iP;out.width(16);out<<" o x x\n"; iP++; } if (cvN(5)>0){ out.width(1);out<<iP;out.width(16);out<<" x o o\n"; iP++; } if (cvN(6)>0){ out.width(1);out<<iP;out.width(16);out<<" x o x\n"; iP++; } if (cvN(7)>0){ out.width(1);out<<iP;out.width(16);out<<" x x o\n"; iP++; } if (cvN(8)>0){ out.width(1);out<<iP;out.width(16);out<<" x x x\n"; } out<<"o: observed; x: missing\n\n"; if (cvN(1)==0){ cout <<"Sorry, you asked too much for missing data analysis! There is no single case with complete data!"<<endl; out <<"Sorry, you asked too much for missing data analysis!\n There is no single case with complete data!\n"<<endl; exit(1); } if (cvN(8)>0){ cout <<"What do you expect if no data observed for all three variables?"<<endl; out <<"What do you expect if no data observed for all three variables?\n"<<endl; exit(1); } std::string strBoot; cout <<"\nPlease input the bootstrap sample size:"<<endl; out <<"The bootstrap sample size is: "; std::cin>>strBoot; int iBoot=atoi(strBoot.c_str()); out << strBoot<<"\n"; //Running time clock_t tStTime, tEndTime; tStTime=clock(); //cout<<matPatData; int iSize=int(cvN(1)); Matrix mCov=mfCov(matPatData.submatrix(1,iSize,1,3),iSize,3); Matrix mBootStat=mBootCI(matPatData,mCov,cvN,r1,iBoot,alpha); //end time tEndTime=clock(); //Output the missing data pattern out <<"\n\n"; out <<"--------------------------------------------------------"<<"\n"; out <<"| Estimated parameters and confidence interval |"<<"\n"; out <<"--------------------------------------------------------"<<"\n"; out <<"\n"; out.width(35); out<<"Quantile"; out.width(18); out<<"BC"; out.width(20); out<<"BCa";out<<"\n"; out.width(10); out<<"Parameter"; out.width(10);out<<"Estimate"; out.width(8);out<<"L";out.width(10);out<<"U"; out.width(10);out<<"L"; out.width(10);out<<"U";out.width(10);out<<"L";out.width(10);out<<"U";out<<"\n"; int i; out << setiosflags(ios::fixed | ios::showpoint)<< setprecision(4); out.width(10); out << "iM"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(1,i); } out<<"\n"; // output results for intercept of Y out.width(10); out << "iY"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(2,i); } out<<"\n"; // output results for a out.width(10); out << "a"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(3,i); } out<<"\n"; // output results for b out.width(10); out << "b"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(4,i); } out<<"\n"; // output results for c out.width(10); out << "c"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(5,i); } out<<"\n"; // output results for eM2 out.width(10); out << "eM2"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(6,i); } out<<"\n"; // output results for residual errors of Y out.width(10); out << "eY2"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(7,i); } out<<"\n"; // output results for sig1 out.width(10); out << "a*b"; for (i=1; i<=7; i++){ out.width(10); out<<mBootStat(8,i); } //timer out <<"\n\nThe total running time is "<<(double)(tEndTime-tStTime)/CLOCKS_PER_SEC<<" seconds. \n"; cout<<"\nThe parameter estimates and the confidence intervals are:"<<endl; cout<<"Please see the output file [["<<output<<"]] for complete results."<<endl; cout << setiosflags(ios::fixed | ios::showpoint)<< setprecision(4); cout<<setw(10)<<mBootStat; out.close(); return 0; }
lab/v1_on_nov_23_2008.txt · Last modified: 2016/01/24 09:48 by 127.0.0.1