### Site Tools

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 <<"---------------------------------------------------------"<<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 <<"*********************************************************"<<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 <<"----------------------------------------------------------\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 <<"**********************************************************\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

//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;

}