R/bmem.R

In johnnyzhz/bmem: Mediation Analysis with Missing Data Using Bootstrap

# package depended: lavaan, Amelia, parallel, snowfall

bmem.moments<-function(x, type=0){ ##rest 
  #0: listwise deletion
  #1: pairwise deletion
  x<-cbind(1, x)
  n<-nrow(x)
  p<-ncol(x)
  res<-matrix(NA, p, p)	
  x<-as.matrix(x)
  if (type ==0){
    x<-na.omit(x)
    res<-crossprod(x,x)/n 
  }else{
    for (i in 1:p) res[i,i]<-mean( (x[complete.cases(x[,i]),i])^2 )
    for (i in 1:(p-1)){
      for (j in (i+1):p){
        y<-x[,c(i,j)]
        y<-na.omit(y)
        res[i,j]<-res[j,i]<-mean(y[,1]*y[,2])
      }
    }
  }
  rownames(res)<-colnames(res)<-c('(i)', (colnames(x))[2:p]) 
  res[1,1]<-1
  res
}

bmem.list.cov<-function(x){
  ## x: data
  ## moment=FALSE: covariance analysis only
  listdata<-na.omit(x)	
  temp.cov<-cov(listdata)
  temp.cov
}


bmem.pair.cov<-function(x){
  ## x: data
  temp.cov<-cov(x, use='pairwise.complete.obs')
  temp.cov
}

bmem.sem<-function(x, ram, N, ...){
  ##x:cov matrix, ram:model
  
  sem.res<-lavaan::sem(model = ram,sample.cov = x,sample.nobs = N, ...)
  
  lavpartable <- lavaan::partable(sem.res)
  
  model.est<-NA
  model.est<-lavaan::parameterEstimates(sem.res) ##a dataframe, columns:
  ##lhs,op,rhs,label,est,se,z,pvalue,ci.lower,ci.upper
  label <- model.est$label               ##name parameters
  lhs <- model.est$lhs
  op <- model.est$op
  rhs <- model.est$rhs
  for (i in 1:length(label)) {
    if(label[i]==""){
      label[i] <- paste(lhs[i],op[i],rhs[i],collapse = " ")
    }
  }
  est <- model.est$est                  ##est is a vector of evaluation of parameters
  names(est) <- label
  model.sta<-lavaan::fitMeasures(sem.res)       ##model.sta a vector of statistics
  model.fit<-c(model.sta["chisq"],model.sta["gfi"], model.sta["agfi"], model.sta["rmsea"], model.sta["nfi"], model.sta["nnfi"], model.sta["cfi"], model.sta["bic"], model.sta["srmr"])
  
  list(est=est,model.fit=model.fit,lavpartable=lavpartable) 
}

bmem.sobel.ind<-function(sem.object, ind){ ##rest
  est<-sem.object$coeff
  temp<-gsub(' +','',ind)	
  temp<-gsub('-','+', temp, fixed=TRUE)
  temp<-unlist(strsplit(temp, '+', fixed=TRUE))
  m<-length(temp)
  temp.est<-0
  par<-NULL
  for (j in 1:m){
    temp1<-unlist(strsplit(temp[j], '*', fixed=TRUE))
    par<-c(par, temp1)
  }
  ind.exp<-parse(text=ind)
  par.list<-as.list(est[par])	
  est.indirect<-eval(ind.exp, par.list)
  ind.deriv<-deriv(ind.exp, par)
  first.deriv<-eval(ind.deriv, par.list)
  first.deriv<-attributes(first.deriv)$gradient

  if (is.null(sem.object$cov)) {
    var.par<-sem.object$vcov[par,par]
  }else{
    var.par<-sem.object$cov[par,par]
  }
  var.ind<-first.deriv%*%var.par%*%t(first.deriv)
  s.e.ind<-sqrt(var.ind)

  res<-matrix(c(est.indirect, s.e.ind, est.indirect/s.e.ind), 1, 3)
  colnames(res)<-c('Estimate', 'S.E.', 'z-score')
  rownames(res)<-ind
  res
}


bmem.sobel<-function(x, ram,  ...){
  N<-nrow(x)
  temp.cov<-bmem.list.cov(x)
  sem.object<-lavaan::sem(model = ram,sample.cov = temp.cov,sample.nobs = N,...)
  
  all.res <- summary(sem.object)$PE
  
  invisible(list(estimates=all.res))
}

bmem.list<-function(x, ram, ...){
  ##x:data, ram:model
  N<-dim(x)[1]
  temp.cov<-bmem.list.cov(x)
  bmem.sem(temp.cov, ram, N, ...)
}

bmem.pair<-function(x, ram, ...){
  N<-dim(x)[1]
  temp.cov<-bmem.pair.cov(x)
  bmem.sem(temp.cov, ram, N, ...)
}

bmem.mi.cov<-function(x, m=10){
  mi.data<-amelia(x,m,0)$imputations
  mi.cov<-list()
  for (j in 1:m){
    temp.cov<-bmem.pair.cov(mi.data[[j]])
    mi.cov[[j]]<-temp.cov
  }
  mi.cov	
}

bmem.v<-function(x, v){
  x[v,v]
}

bmem.mi<-function(x, ram, v, m=10, ...){
  N<-dim(x)[1]
  if (missing(v)) v<-1:(ncol(x))
  temp.cov<-bmem.mi.cov(x,m)
  mi.est<-NULL
  mi.fit<-NULL
  for (j in 1:m){
    v.cov<-bmem.v(temp.cov[[j]], v)
    temp<-bmem.sem(v.cov, ram, N, ...)
    mi.est<-rbind(mi.est, temp$est)
    mi.fit<-rbind(mi.fit, temp$model.fit)
  }
  mi.est<-apply(mi.est,2,mean)
  mi.fit<-apply(mi.fit,2,mean)
  list(est=mi.est, model.fit=mi.fit,lavpartable=temp$lavpartable)
}

bmem.pattern<-function(x){
  if (!is.matrix(x)) x<-data.matrix(x)
  n<-dim(x)[1]
  p<-dim(x)[2]
  misorder<-rep(0,n)
  for (i in 1:n){
    misorderj<-0
    for (j in 1:p){
      if (is.na(x[i,j])) misorderj<-misorderj+2^(j-1)
    }
    misorder[i]<-misorderj
  }
  temp<-order(misorder)
  x<-x[temp,]
  misn<-misorder[temp]
  
  mi<-0; nmi<-0;oi<-0; noi<-0;
  for (j in 1:p){
    if (is.na(x[1,j])){
      mi<-c(mi,j)  ## recording the missing variable subscript in the first case
      nmi<-nmi+1   ## number of missing values in the first case
    }else{
      oi<-c(oi,j)
      noi<-noi+1
    }
  }
  oi<-oi[2:(noi+1)]
  if (nmi==0){
    misinfo_0 = c(noi, oi)
  }else{
    mi<-mi[2:(nmi+1)]
    misinfo_0<-c(noi,oi,mi) ##recording the # observed variables, locations;
  }	
  patnobs <- 0 ## number of observed cases in a pattern
  totpat<-1; ncount<-1;
  t1<-misn[1]
  for (i in 2:n){
    if (misn[i]==t1){
      ncount<-ncount+1
    }else{
      patnobs<-c(patnobs,ncount)
      t1<-misn[i]
      ncount<-1
      totpat<-totpat+1
      mi<-0; nmi<-0;oi<-0; noi<-0;
      for (j in 1:p){
        if (is.na(x[i,j])){
          mi<-c(mi,j)
          nmi<-nmi+1
        }else{
          oi<-c(oi,j)
          noi<-noi+1
        }
      }
      oi<-oi[2:(noi+1)]
      mi<-mi[2:(nmi+1)]
      misinfo_0 <- rbind(misinfo_0, c(noi,oi,mi))
    }
  }
  patnobs<-c(patnobs, ncount)
  patnobs<-patnobs[2:(totpat+1)]
  if (is.vector(misinfo_0)) {misinfo<-c(patnobs, misinfo_0)}else{misinfo<-cbind(patnobs, misinfo_0)}
  if (!is.matrix(misinfo)){misinfo<-matrix(misinfo, nrow=1)}
  colnames(misinfo)<-NULL
  rownames(misinfo)<-NULL
  list(misinfo=misinfo, x=x)
}

bmem.ssq<-function(x){
  sum(x^2)
}

bmem.em.cov<-function(xmis, max_it=500){
  x<-xmis[[2]]
  misinfo<-xmis[[1]]
  ep <- 1e-12  ## precision
  n<-dim(x)[1]
  p<-dim(x)[2]
  mu0<-apply(x,2,mean, na.rm=TRUE) ## starting values for mu
  sig0<-cov(x,use="complete.obs") ## starting values for sigma
  n_it<-0;        
  err<-0
  dt<-1;
  while (dt>ep && n_it <= max_it){
    sumx<-rep(0,p); sumxx<-array(0,dim=c(p,p)); sumw1<-0; sumw2<-0;
    npat<-dim(misinfo)[1]  ## number of missing data patterns
    p1<-misinfo[1,2]       ## number of observed variables in pattern 1
    n1<-misinfo[1,1]       ## number of cases in pattern 1
    if (p1==p){            ## complete data
      sigin <- solve(sig0)  ## matrix inverse
      for (i in 1:n1){
        xi<-x[i,]
        xi0<-xi-mu0
        sumw1<-sumw1+1;
        xxi0<-xi0%*%t(xi0)
        sumx<-sumx+xi
        sumxx<-sumxx+xxi0
        sumw2<-sumw2+1
      } ## end for
    }else{ ## end p1==p
      ## with missing data			
      o1<-misinfo[1,3:(2+p1)]
      m1<-misinfo[1,(2+p1+1):(p+2)]
      mu_o<-mu0[o1]; mu_m<-mu0[m1]
      sig_oo<-sig0[o1,o1]; sig_om<-sig0[o1,m1];
      if (p1==1) {sig_mo<-sig_om}else{sig_mo<-t(sig_om)}
      sig_mm<-sig0[m1,m1];
      sigin_oo<-solve(sig_oo)
      beta_mo<-sig_mo%*%sigin_oo
      
      delt <- array(0, dim=c(p,p))
      delt[m1,m1]<-sig_mm - beta_mo%*%sig_om
      for (i in 1:n1){
        xi<-x[i,]
        xi_o<-xi[o1]
        xi0_o<-xi_o-mu_o
        stdxi_o<-sigin_oo%*%xi0_o
        sumw1<-sumw1+1
        xm1<-mu_m+sig_mo%*%stdxi_o
        xi[m1]<-xm1
        xi0<-xi-mu0
        xxi0<-xi0%*%t(xi0)
        sumx<-sumx+xi
        sumxx<-sumxx+xxi0+delt
        sumw2<-sumw2+1
      } ##end for 1:n1  
    }## end of (p1=p)
    ## start from pattern 2	
    if (npat>1){
      snj<-n1	
      for (j in 2:npat){
        nj<-misinfo[j,1]; pj<-misinfo[j,2];
        oj<-misinfo[j, 3:(2+pj)]; mj<-misinfo[j, (2+pj+1):(p+2)];
        mu_o<-mu0[oj]; mu_m<-mu0[mj];
        sig_oo<-sig0[oj,oj]; sig_om<-sig0[oj,mj];
        if (pj==1) {sig_mo<-sig_om}else{sig_mo<-t(sig_om)}  
        sig_mm<-sig0[mj,mj];
        sigin_oo<-solve(sig_oo)
        beta_mo<-sig_mo%*%sigin_oo
        delt <- array(0, dim=c(p,p))
        delt[mj,mj]<-sig_mm - beta_mo%*%sig_om
        
        for (i in ((snj+1):(snj+nj))){
          xi<-x[i,]
          xi_o<-xi[oj]
          xi0_o<-xi_o - mu_o
          stdxi_o<-sigin_oo%*%xi0_o					
          sumw1<-sumw1+1
          xmj<-mu_m+sig_mo%*%stdxi_o
          xi[mj]<-xmj
          xi0<-xi-mu0
          xxi0<-xi0%*%t(xi0)
          sumx<-sumx+1*xi
          sumxx<-sumxx+xxi0+delt
          sumw2<-sumw2+1
        }
        snj<-snj+nj
      } ## for (j in 2:npat)
    }
    mu1<-sumx/sumw1
    sig1<-sumxx/n
    dt<-(bmem.ssq(mu1-mu0)+bmem.ssq(sig1-sig0))/(bmem.ssq(mu0)+bmem.ssq(sig0));
    mu0<-mu1;
    sig0<-sig1;
    n_it<-n_it+1;
  } ## end while
  if (n_it == max_it) warnings('Maximum iteration for EM algorithm is exceeded!')
  temp.cov<-sig1
  rownames(temp.cov)<-colnames(sig1)
  m<-matrix(mu1, p, 1)
  temp.cov
}

bmem.em.rcov<-function(xmis, varphi=.1, max_it=1000, st='i'){
  ## x: data set
  ## misinfo: missing data pattern
  ## varphi: 
  x<-xmis[[2]]
  misinfo<-xmis[[1]]
  ep <- 1e-6  ## precision
  n<-dim(x)[1]
  p<-dim(x)[2]
  
  mu0<-rep(0,p)
  sig0<-diag(p)
  if (st=='mcd'){
    y<-na.omit(x)
    ny<-nrow(y)
    par.st<-cov.rob(y, method='mcd')
    mu0<-par.st$center
    sig0<-par.st$cov
  }
  
  n_it<-0;        
  dt<-1;
  if (varphi==0){
    ck<-10e+10
    cbeta<-1
  }else{
    prob<-1-varphi ## chi-square p-value
    chip<-qchisq(prob, p)
    ck<-sqrt(chip)
    cbeta<-( p*pchisq(chip, p+2) + chip*(1-prob) )/p
  }
  while (dt>ep && n_it <= max_it){
    sumx<-rep(0,p); sumxx<-array(0,dim=c(p,p)); sumw1<-0; sumw2<-0;
    npat<-dim(misinfo)[1]  ## number of missing data patterns
    p1<-misinfo[1,2]       ## number of observed variables in pattern 1
    n1<-misinfo[1,1]       ## number of cases in pattern 1
    if (p1==p){            ## complete data
      sigin <- solve(sig0)  ## matrix inverse
      for (i in 1:n1){
        xi<-x[i,]
        xi0<-xi-mu0
        di2<-xi0%*%sigin%*%xi0
        di<-sqrt(di2)
        ## Huber weight functions
        if (di<=ck){
          wi1<-1
          wi2<-1/cbeta
        }else{
          wi1<-ck/di
          wi2<-wi1*wi1/cbeta
        } ## end Huber weight
        sumw1<-sumw1+wi1;
        xxi0<-xi0%*%t(xi0)
        sumx<-sumx+wi1*xi
        sumxx<-sumxx+c(wi2)*xxi0
        sumw2<-sumw2+wi2
      } ## end for
    }else{ ## end p1==p
      ## with missing data
      if (varphi==0){
        ck1<-1e+10
        cbeta1<-1
      }else{ 
        chip1<-qchisq(prob, p1)
        ck1<-sqrt(chip1)
        cbeta1<-( p1*pchisq(chip1,p1+2) + chip1*(1-prob) )/p1
      }
      o1<-misinfo[1,3:(2+p1)]
      m1<-misinfo[1,(2+p1+1):(p+2)]
      mu_o<-mu0[o1]; mu_m<-mu0[m1]
      sig_oo<-sig0[o1,o1]; sig_om<-sig0[o1,m1];
      if (p1==1) {sig_mo<-sig_om}else{sig_mo<-t(sig_om)} 
      sig_mm<-sig0[m1,m1];
      sigin_oo<-solve(sig_oo)
      beta_mo<-sig_mo%*%sigin_oo
      
      delt <- array(0, dim=c(p,p))
      delt[m1,m1]<-sig_mm - beta_mo%*%sig_om
      for (i in 1:n1){
        xi<-x[i,]
        xi_o<-xi[o1]
        xi0_o<-xi_o-mu_o
        stdxi_o<-sigin_oo%*%xi0_o
        di2<-xi0_o%*%stdxi_o
        di<-sqrt(di2)
        if (di<=ck1){ ##Huber weight
          wi1<-1
          wi2<-1/cbeta1
        }else{
          wi1<-ck1/di
          wi2<-wi1*wi1/cbeta1
        }
        sumw1<-sumw1+wi1
        xm1<-mu_m+sig_mo%*%stdxi_o
        xi[m1]<-xm1
        xi0<-xi-mu0
        xxi0<-xi0%*%t(xi0)
        sumx<-sumx+wi1*xi
        sumxx<-sumxx+c(wi2)*xxi0+delt
        sumw2<-sumw2+wi2
      } ##end for 1:n1  
    }## end of (p1=p)
    ## start from pattern 2	
    if (npat>1){
      snj<-n1	
      for (j in 2:npat){
        nj<-misinfo[j,1]; pj<-misinfo[j,2];
        oj<-misinfo[j, 3:(2+pj)]; mj<-misinfo[j, (2+pj+1):(p+2)];
        mu_o<-mu0[oj]; mu_m<-mu0[mj];
        sig_oo<-sig0[oj,oj]; sig_om<-sig0[oj,mj];
        if (pj==1) {sig_mo<-sig_om}else{sig_mo<-t(sig_om)} 
        sig_mm<-sig0[mj,mj];
        sigin_oo<-solve(sig_oo)
        beta_mo<-sig_mo%*%sigin_oo
        delt <- array(0, dim=c(p,p))
        delt[mj,mj]<-sig_mm - beta_mo%*%sig_om
        if (varphi==0){
          ckj<-10e+10
          cbetaj<-1
        }else{
          chipj<-qchisq(prob,pj)
          ckj<-sqrt(chipj)
          cbetaj<- ( pj*pchisq(chipj, pj+2) + chipj*(1-prob) )/pj
        }
        for (i in ((snj+1):(snj+nj))){
          xi<-x[i,]
          xi_o<-xi[oj]
          xi0_o<-xi_o - mu_o
          stdxi_o<-sigin_oo%*%xi0_o
          di2<-xi0_o%*%stdxi_o
          di<-sqrt(di2)
          if (di<=ckj){ ##Huber weight
            wi1<-1
            wi2<-1/cbetaj
          }else{
            wi1<-ckj/di
            wi2<-wi1*wi1/cbetaj
          }
          sumw1<-sumw1+wi1
          xmj<-mu_m+sig_mo%*%stdxi_o
          xi[mj]<-xmj
          xi0<-xi-mu0
          xxi0<-xi0%*%t(xi0)
          sumx<-sumx+wi1*xi
          sumxx<-sumxx+c(wi2)*xxi0+delt
          sumw2<-sumw2+wi2
        }
        snj<-snj+nj
      } ## for (j in 2:npat)
    }
    mu1<-sumx/sumw1
    sig1<-sumxx/n
    dt<-max(c(max(abs(mu1-mu0)), max(abs(sig1-sig0))));
    mu0<-mu1;
    sig0<-sig1;
    n_it<-n_it+1;
  } ## end while
  if (n_it == max_it) warnings('Maximum iteration for EM algorithm is exceeded!')
  temp.cov<-sig1
  rownames(temp.cov)<-colnames(temp.cov)<-colnames(x)
  m<-matrix(mu1, p, 1)
  temp.cov
}

bmem.em<-function(x, ram, v, robust=FALSE, varphi=.1, st='i', max_it=500, ...){
  N<-dim(x)[1]
  if (missing(v)) v<-1:(ncol(x))
  temp.pattern<-bmem.pattern(x)
  if (robust){
    temp.cov<-bmem.em.rcov(temp.pattern, varphi, max_it, st)
  }else{
    temp.cov<-bmem.em.cov(temp.pattern, max_it)
  }
  temp.cov<-bmem.v(temp.cov,v)
  bmem.sem(temp.cov, ram, N, ...)
}

bmem.list.jack<-function(x, ram, ...){
  n<-dim(x)[1]
  nseq<-1:n
  jack.est<-NULL
  jack.fit<-NULL
  for (i in 1:n){
    x.jack<-x[nseq[-i],]
    jack.temp<-try(bmem.list(x.jack, ram, ...))
    if (class(jack.temp)!="try-error"){
      jack.est<-rbind(jack.est, jack.temp$est)
      jack.fit<-rbind(jack.fit, jack.temp$model.fit)
    }
  }
  list(jack.est=jack.est, jack.fit=jack.fit,lavpartable=jack.temp$lavpartable)
}

bmem.pair.jack<-function(x, ram, ...){
  n<-dim(x)[1]
  nseq<-1:n
  jack.est<-NULL
  jack.fit<-NULL
  for (i in 1:n){
    x.jack<-x[nseq[-i],]
    jack.temp<-try(bmem.pair(x.jack, ram, ...))
    if (class(jack.temp)!="try-error"){
      jack.est<-rbind(jack.est, jack.temp$est)
      jack.fit<-rbind(jack.fit, jack.temp$model.fit)
    }
  }
  list(jack.est=jack.est, jack.fit=jack.fit,lavpartable=jack.temp$lavpartable)
}


bmem.mi.jack<-function(x, ram, v, m=10, ...){
  n<-dim(x)[1]
  nseq<-1:n
  jack.est<-NULL
  jack.fit<-NULL
  for (i in 1:n){
    x.jack<-x[nseq[-i],]
    jack.temp<-try(bmem.mi(x.jack, ram, v, m, ...))
    if (class(jack.temp)!="try-error"){
      jack.est<-rbind(jack.est, jack.temp$est)
      jack.fit<-rbind(jack.fit, jack.temp$model.fit)
    }
  }
  list(jack.est=jack.est, jack.fit=jack.fit,lavpartable=jack.temp$lavpartable)
}

bmem.em.jack<-function(x, ram, v, robust=FALSE, varphi=.1, st='i', max_it=500, ...){  
  n<-dim(x)[1]
  nseq<-1:n
  jack.est<-NULL
  jack.fit<-NULL
  for (i in 1:n){
    x.jack<-x[nseq[-i],]
    jack.temp<-try(bmem.em(x.jack, ram, v, robust, varphi, st, max_it, ...))
    if (class(jack.temp)!="try-error"){
      jack.est<-rbind(jack.est, jack.temp$est)
      jack.fit<-rbind(jack.fit, jack.temp$model.fit)
    }
  }
  list(jack.est=jack.est, jack.fit=jack.fit,lavpartable=jack.temp$lavpartable)
}

bmem.list.boot<-function(x, ram, boot=1000, parallel=FALSE,ncore = 1,...){	
  # x:data, ram:model
  model0<-bmem.list(x, ram, ...)
  par0<-model0$est
  fit0<-model0$model.fit
  n<-dim(x)[1]
  
  myfun <- function(i){
    x.boot<-x[sample(n,n, replace=TRUE),]
    options(warn =-1)  ## hide all warnings
    modelb<-try(bmem.list(x.boot, ram, ...),silent=TRUE) ## silent=TRUE
    if (class(modelb)!="try-error"){
      return(list(est=modelb$est,fit=modelb$model.fit))  ## use result[[i]]$est,result[[i]]$fit to get the ith result  
    }
  }
  old_options <- options(); options(warn = -1)
  if(parallel){ ## needs package "snowfall" and "parallel"
    snowfall::sfInit(parallel = TRUE, cpus = ncore)# detectCores() - 1)
    snowfall::sfLibrary("lavaan", character.only = TRUE)
    sfExport("bmem.list","bmem.list.cov","bmem.sem")
    
    result <- snowfall::sfLapply(1:boot, myfun) 
    sfStop()
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:length(result)) {
      boot.est<-rbind(boot.est, result[[i]]$est)
      boot.fit<-rbind(boot.fit, result[[i]]$fit)
    }
  }else{
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:boot){
      x.boot<-x[sample(n,n, replace=TRUE),]
      options(warn =-1)  ##hide all warnings
      modelb<-try(bmem.list(x.boot, ram, ...),silent=TRUE) ##silent=TRUE
      if (class(modelb)!="try-error"){
        boot.est<-rbind(boot.est, modelb$est)
        boot.fit<-rbind(boot.fit, modelb$model.fit)
      }
    }
  }
  options(old_options)
  colnames(boot.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')	

  list(par.boot=boot.est, par0=par0, boot.fit=boot.fit, fit0=fit0, lavpartable=model0$lavpartable) ##"0" is the result of original data

}

bmem.pair.boot<-function(x, ram, boot=1000, parallel=FALSE,ncore = 1,...){
  model0<-bmem.pair(x, ram, ...)
  par0<-model0$est
  fit0<-model0$model.fit
  n<-dim(x)[1]
  
  myfun <- function(i){
    x.boot<-x[sample(n,n, replace=TRUE),]
    options(warn =-1)  ## hide all warnings
    modelb<-try(bmem.pair(x.boot, ram, ...),silent=TRUE) ## silent=TRUE
    if (class(modelb)!="try-error"){
      return(list(est=modelb$est,fit=modelb$model.fit))  ## use result[[i]]$est,result[[i]]$fit to get the ith result  
    }
  }
  old_options <- options(); options(warn = -1)
  if(parallel){ ## needs package "snowfall" and "parallel"
    snowfall::sfInit(parallel = TRUE, cpus = ncore)# detectCores() - 1)
    snowfall::sfLibrary("lavaan", character.only = TRUE)
    sfExport("bmem.pair","bmem.pair.cov","bmem.sem")
    
    result <- snowfall::sfLapply(1:boot, myfun) 
    sfStop()
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:length(result)) {
      boot.est<-rbind(boot.est, result[[i]]$est)
      boot.fit<-rbind(boot.fit, result[[i]]$fit)
    }
  }else{
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:boot){
      x.boot<-x[sample(n,n, replace=TRUE),]
      options(warn =-1)  ##hide all warnings
      modelb<-try(bmem.pair(x.boot, ram, ...),silent=TRUE) ##silent=TRUE
      if (class(modelb)!="try-error"){
        boot.est<-rbind(boot.est, modelb$est)
        boot.fit<-rbind(boot.fit, modelb$model.fit)
      }
    }
  }
  options(old_options)
  
  colnames(boot.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')	
  list(par.boot=boot.est, par0=par0, boot.fit=boot.fit, fit0=fit0,lavpartable=model0$lavpartable)
}

bmem.mi.boot<-function(x, ram, v, m=10, boot=1000, parallel=FALSE,ncore = 1,...){
  model0<-bmem.mi(x, ram, v, m, ...)
  par0<-model0$est
  fit0<-model0$model.fit
  n<-dim(x)[1]
  
  myfun <- function(i){
    x.boot<-x[sample(n,n, replace=TRUE),]
    options(warn =-1)  ## hide all warnings
    modelb<-try(bmem.mi(x.boot, ram, ...),silent=TRUE) ## silent=TRUE
    if (class(modelb)!="try-error"){
      return(list(est=modelb$est,fit=modelb$model.fit))  ## use result[[i]]$est,result[[i]]$fit to get the ith result  
    }
  }
  old_options <- options(); options(warn = -1)
  if(parallel){ ## needs package "snowfall" and "parallel"
    snowfall::sfInit(parallel = TRUE, cpus = ncore)# detectCores() - 1)
    snowfall::sfLibrary("lavaan", character.only = TRUE)
    sfExport("bmem.mi","bmem.mi.cov","bmem.sem")
    
    result <- snowfall::sfLapply(1:boot, myfun) 
    sfStop()
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:length(result)) {
      boot.est<-rbind(boot.est, result[[i]]$est)
      boot.fit<-rbind(boot.fit, result[[i]]$fit)
    }
  }else{
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:boot){
      x.boot<-x[sample(n,n, replace=TRUE),]
      options(warn =-1)  ##hide all warnings
      modelb<-try(bmem.mi(x.boot, ram, ...),silent=TRUE) ##silent=TRUE
      if (class(modelb)!="try-error"){
        boot.est<-rbind(boot.est, modelb$est)
        boot.fit<-rbind(boot.fit, modelb$model.fit)
      }
    }
  }
  options(old_options)
  
  colnames(boot.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')		
  list(par.boot=boot.est, par0=par0, boot.fit=boot.fit, fit0=fit0, lavpartable=model0$lavpartable)
}


bmem.em.boot<-function(x, ram, v, robust=FALSE, varphi=.1, st='i', boot=1000, max_it=500, parallel=FALSE,ncore = 1,...){
  model0<-bmem.em(x, ram, v, robust, varphi, st, max_it, ...)
  par0<-model0$est
  fit0<-model0$model.fit
  n<-dim(x)[1]
  
  myfun <- function(i){
    x.boot<-x[sample(n,n, replace=TRUE),]
    options(warn =-1)  ## hide all warnings
    modelb<-try(bmem.em(x.boot, ram, ...),silent=TRUE) ## silent=TRUE
    if (class(modelb)!="try-error"){
      return(list(est=modelb$est,fit=modelb$model.fit))  ## use result[[i]]$est,result[[i]]$fit to get the ith result  
    }
  }
  old_options <- options(); options(warn = -1)
  if(parallel){ ## needs package "snowfall" and "parallel"
    snowfall::sfInit(parallel = TRUE, cpus = ncore)# detectCores() - 1)
    snowfall::sfLibrary("lavaan", character.only = TRUE)
    sfExport("bmem.em","bmem.em.cov","bmem.sem")
    
    result <- snowfall::sfLapply(1:boot, myfun) 
    sfStop()
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:length(result)) {
      boot.est<-rbind(boot.est, result[[i]]$est)
      boot.fit<-rbind(boot.fit, result[[i]]$fit)
    }
  }else{
    
    boot.est<-NULL
    boot.fit<-NULL
    for (i in 1:boot){
      x.boot<-x[sample(n,n, replace=TRUE),]
      options(warn =-1)  ##hide all warnings
      modelb<-try(bmem.em(x.boot, ram, ...),silent=TRUE) ##silent=TRUE
      if (class(modelb)!="try-error"){
        boot.est<-rbind(boot.est, modelb$est)
        boot.fit<-rbind(boot.fit, modelb$model.fit)
      }
    }
  }
  options(old_options)
  
  colnames(boot.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')	
  list(par.boot=boot.est, par0=par0, boot.fit=boot.fit, fit0=fit0, lavpartable=model0$lavpartable)
}

## bootstrap confidence intervals
bmem.ci.p<-function(par.boot, par0, cl=.95){
  alpha<-(1-cl)/2
  alpha<-c(alpha, 1-alpha)
  alpha<-sort(alpha)
  ci<-apply(par.boot, 2, quantile, prob=alpha, na.rm=TRUE)
  se.boot<-apply(par.boot, 2, sd, na.rm=TRUE)
  estimate<-par0
  cbind(estimate, se.boot, t(ci))
}

bmem.ci.norm<-function(par.boot, par0, cl=.95){
  alpha<-(1-cl)/2
  alpha<-c(alpha, 1-alpha)
  alpha<-sort(alpha)
  se.boot<-apply(par.boot, 2, sd, na.rm=TRUE)
  estimate<-par0
  ci<-estimate+qnorm(alpha)%*%t(se.boot)
  dig <- max(2L, getOption("digits"))
  np<-length(alpha)
  qs <- paste(if (np < 100) 
    formatC(100 * alpha, format = "fg", width = 1, digits = dig)
    else format(100 * alpha, trim = TRUE, digits = dig), 
    "%", sep = "")
  ci.res<-cbind(estimate, se.boot, t(ci))
  colnames(ci.res)<-c('estimate','se.boot', qs)
  ci.res
}

bmem.ci.bc1<-function(x, b, cl=.95){
  n<-length(x)
  z0<-qnorm(sum(x<b, na.rm=TRUE)/n)
  alpha<-(1-cl)/2
  alpha<-c(alpha, 1-alpha)
  alpha<-sort(alpha)
  alpha1<-alpha
  alpha<-pnorm(2*z0+qnorm(alpha))
  dig <- max(2L, getOption("digits"))
  np<-length(alpha)
  qs<-quantile(x, alpha, na.rm=TRUE)
  names(qs) <- paste(if (np < 100) 
    formatC(100 * alpha1, format = "fg", width = 1, digits = dig)
    else format(100 * alpha1, trim = TRUE, digits = dig), 
    "%", sep = "")
  qs
}

bmem.ci.bc<-function(par.boot, par0, cl=.95){
  se.boot<-apply(par.boot, 2, sd, na.rm=TRUE)
  estimate<-par0
  p<-ncol(par.boot)
  ci<-NULL
  for (i in 1:p){
    ci<-rbind(ci, bmem.ci.bc1(par.boot[,i], par0[i], cl))
  }
  cbind(estimate, se.boot, ci)
}

bmem.ci.bca1<-function(x, b, jack, cl=.95){
  n<-length(x)
  z0<-qnorm(sum(x<b, na.rm=TRUE)/n)
  alpha<-(1-cl)/2
  alpha<-c(alpha, 1-alpha)
  alpha<-sort(alpha)
  alpha1<-alpha
  mjack<-mean(jack, na.rm=TRUE)
  a<-sum((mjack-jack)^3, na.rm=TRUE)/(6*(sum((mjack-jack)^2, na.rm=TRUE))^1.5)
  alpha<-pnorm(z0+(z0+qnorm(alpha))/(1-a*(z0+qnorm(alpha))))
  dig <- max(2L, getOption("digits"))
  np<-length(alpha)
  if (alpha[1]=='NaN'){
    qs<-c(NA,NA)
  }else{
    qs<-quantile(x, alpha, na.rm=TRUE)
  }
  names(qs) <- paste(if (np < 100) 
    formatC(100 * alpha1, format = "fg", width = 1, digits = dig)
    else format(100 * alpha1, trim = TRUE, digits = dig), 
    "%", sep = "")
  qs
}

bmem.ci.bca<-function(par.boot, par0, jack, cl=.95){
  se.boot<-apply(par.boot, 2, sd)
  estimate<-par0
  p<-ncol(par.boot)
  ci<-NULL
  for (i in 1:p){
    ci<-rbind(ci, bmem.ci.bca1(par.boot[,i], par0[i], jack[,i], cl))
  }
  cbind(estimate, se.boot, ci)
}


### A main function for analysis
bmem<-function(data, model, v, method='tsml', ci='bc', cl=.95, boot=1000, m=10, varphi=.1, st='i', robust=FALSE, max_it=500, parallel=FALSE, ncore=1, ...){
  ## data:data, model:model
  ## method: list, pair, mi, tsml
  ## ci: norm, perc, bc, bca
  N<-nrow(data)
  P<-ncol(data)
  if (missing(v)) v<-1:P
  
  ## for listwise deletion method
  if (method=='list'){
    data<-data[,v]
    boot.est<-bmem.list.boot(data, model, boot, parallel, ncore, ...)
    if (ci=='norm'){
      ci.est<-bmem.ci.norm(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.norm(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='perc'){ 
      ci.est<-bmem.ci.p(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.p(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bc'){ 
      ci.est<-bmem.ci.bc(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.bc(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bca'){
      jack.est<-bmem.list.jack(data,model, ...)
      ci.est<-bmem.ci.bca(boot.est$par.boot, boot.est$par0, jack.est$jack.est, cl)
      ci.fit<-bmem.ci.bca(boot.est$boot.fit, boot.est$fit0, jack.est$jack.fit, cl)
    }
  }
  
  ## for pairwise deletion method
  if (method=='pair'){
    data<-data[,v]
    boot.est<-bmem.pair.boot(data, model, boot, parallel, ncore,...)
    if (ci=='norm'){
      ci.est<-bmem.ci.norm(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.norm(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='perc'){ 
      ci.est<-bmem.ci.p(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.p(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bc'){ 
      ci.est<-bmem.ci.bc(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.bc(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bca'){
      jack.est<-bmem.pair.jack(data,model, ...)
      ci.est<-bmem.ci.bca(boot.est$par.boot, boot.est$par0, jack.est$jack.est, cl)
      ci.fit<-bmem.ci.bca(boot.est$boot.fit, boot.est$fit0, jack.est$jack.fit, cl)
    }
  }
  
  ## for multiple imputation method
  if (method=='mi'){
    boot.est<-bmem.mi.boot(data,model,v,m,boot,parallel, ncore,...)
    if (ci=='norm'){
      ci.est<-bmem.ci.norm(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.norm(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='perc'){ 
      ci.est<-bmem.ci.p(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.p(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bc'){ 
      ci.est<-bmem.ci.bc(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.bc(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bca'){
      jack.est<-bmem.mi.jack(data,model,v,m,...)
      ci.est<-bmem.ci.bca(boot.est$par.boot, boot.est$par0, jack.est$jack.est, cl)
      ci.fit<-bmem.ci.bca(boot.est$boot.fit, boot.est$fit0, jack.est$jack.fit, cl)
    }
  }
  
  ## for EM method
  if (method=='tsml'){
    boot.est<-bmem.em.boot(data,model, v, robust, varphi, st, boot, max_it, parallel, ncore,...)
    if (ci=='norm'){
      ci.est<-bmem.ci.norm(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.norm(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='perc'){ 
      ci.est<-bmem.ci.p(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.p(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bc'){ 
      ci.est<-bmem.ci.bc(boot.est$par.boot, boot.est$par0, cl)
      ci.fit<-bmem.ci.bc(boot.est$boot.fit, boot.est$fit0, cl)
    }
    if (ci=='bca'){
      jack.est<-bmem.em.jack(data,model, v, robust, varphi, st, max_it, ...)
      ci.est<-bmem.ci.bca(boot.est$par.boot, boot.est$par0, jack.est$jack.est, cl)
      ci.fit<-bmem.ci.bca(boot.est$boot.fit, boot.est$fit0, jack.est$jack.fit, cl)
    }
  }
  cat('The bootstrap confidence intervals for parameter estimates\n')
  print(ci.est)
  
  cat('\nThe bootstrap confidence intervals for model fit indices\n')
  rownames(ci.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')
  print(ci.fit)
  cat('\nThe literature has suggested the use of Bollen-Stine bootstrap for model fit. To do so, use the function bmem.bs().\n')
  
  lavpartable <- boot.est$lavpartable[,c(2:4,7,10,11,14,15)]## to be written
  lavpartable[,"est"] <- ci.est[,"estimate"]
  lavpartable[,"se"] <- ci.est[,"se.boot"]
  
  infor <- list(samplesize=N,method=method,boot=boot,sucboot=nrow(boot.est$par.boot),ci=ci)
  
  if (ci=='bca') {
    bmemobject<-list(ci=ci.est, ci.fit=ci.fit, boot.est=boot.est, jack.est=jack.est, lavpartable=lavpartable,infor=infor)
  }else{
    bmemobject<-list(ci=ci.est, ci.fit=ci.fit, boot.est=boot.est, lavpartable=lavpartable,infor=infor)
  }
  class(bmemobject)<-'bmem'
  invisible(bmemobject)
}

bmem.raw2cov<-function(x){
  p<-nrow(x)
  m<-as.matrix(x[1, 2:p])
  Sraw<-x[2:p, 2:p]
  S<-Sraw-m%*%t(m)
  list(S=S, m=m)
}

bmem.bs<-function(x, ram, v, ci='bc', cl=.95, boot=1000, max_it=500, ...){
  ## Estimate the saturated mean and covariance matrix
  xmiss<-bmem.pattern(x)
  s.cov<-bmem.em.cov(xmiss, max_it)
  
  ## Estimate the model indicated covariance matrix
  N<-nrow(x)
  temp.sem<-sem(ram, s.cov, N, raw=TRUE, fixed.x='(i)', ...)
  m.cov<-temp.sem$C
  
  ## Take out the variables in the model (excluding the auxiliary variables)
  obsvar<-rownames(m.cov)
  s.cov.m<-s.cov[obsvar, obsvar]
  
  ## mean and covariance matrix
  s.mcov<-bmem.raw2cov(s.cov.m)
  m.mcov<-bmem.raw2cov(m.cov)
  
  ## select data
  x.m<-xmiss$x[,obsvar[2:length(obsvar)]]
  
  ## generate new data based on the model
  x.bs<-x.m
  x.m.miss<-bmem.pattern(x.m)
  
  npat<-nrow(x.m.miss$misinfo)
  snj<-0
  for (j in 1:npat){
    nj<-x.m.miss$misinfo[j,1]
    pj<-x.m.miss$misinfo[j,2]
    oj<-x.m.miss$misinfo[j, 3:(2+pj)]
    mu_m<-m.mcov$m[oj,1]
    sig_m<-m.mcov$S[oj,oj]
    mu_s<-s.mcov$m[oj,1]
    sig_s<-s.mcov$S[oj,oj]
    for (i in (snj+1):(snj+nj)){
      x.bs[i, oj]<-(x.m.miss$x[i,oj] - t(mu_s)%*%solve(sig_s)%*%sig_m+t(mu_m))
    }
    snj<-nj	
  }
  
  ## combine model data with auxiliary variables
  obsvar<-obsvar[-1]
  x[, obsvar]<-x.bs
  
  ## now for bootstrap based on the newly generated sample
  N<-nrow(x)
  P<-ncol(x)
  if (missing(v)) v<-1:P
  ##boot.est<-bmem.em.boot(x,ram,indirect,v,boot,moment,max_it, ...)
  boot.est<-bmem.em.boot(x,ram,v,boot,max_it, ...)
  if (ci=='norm'){
    ci.est<-bmem.ci.norm(boot.est$par.boot, boot.est$par0, cl)
    ci.fit<-bmem.ci.norm(boot.est$boot.fit, boot.est$fit0, cl)
  }
  if (ci=='perc'){ 
    ci.est<-bmem.ci.p(boot.est$par.boot, boot.est$par0, cl)
    ci.fit<-bmem.ci.p(boot.est$boot.fit, boot.est$fit0, cl)
  }
  if (ci=='bc'){ 
    ci.est<-bmem.ci.bc(boot.est$par.boot, boot.est$par0, cl)
    ci.fit<-bmem.ci.bc(boot.est$boot.fit, boot.est$fit0, cl)
  }
  if (ci=='bca'){
    jack.est<-bmem.em.jack(x,ram,v,max_it, ...)
    ci.est<-bmem.ci.bca(boot.est$par.boot, boot.est$par0, jack.est$jack.est, cl)
    ci.fit<-bmem.ci.bca(boot.est$boot.fit, boot.est$fit0, jack.est$jack.fit, cl)
  }
  cat('The bootstrap confidence intervals for parameter estimates\n')
  print(ci.est)
  
  cat('\nThe bootstrap confidence intervals for model fit indices\n')
  rownames(ci.fit)<-rownames(ci.fit)<-c('chisq', 'GFI','AGFI', 'RMSEA','NFI','NNFI','CFI','BIC','SRMR')
  print(ci.fit)
  cat('\nThe literature has suggested the use of raw data bootstrap for standard errors. To do so, use the function bmem().\n')
  
  if (ci=='bca') {
    bmemobject<-list(ci=ci.est, ci.fit=ci.fit, boot.est=boot.est, jack.est=jack.est)
  }else{
    bmemobject<-list(ci=ci.est, ci.fit=ci.fit, boot.est=boot.est)
  }
  class(bmemobject)<-'bmem'
  invisible(bmemobject)
}

bmem.plot<-function(x, par,...){
  ## x: output from bmem function
  ## par: a parameter or a fit indice to use
  parnames<-colnames(x$boot.est$par.boot)
  fitnames<-colnames(x$boot.est$boot.fit)
  
  if (par %in% parnames){
    hist(x$boot.est$par.boot[,par], xlab=par, ylab='Density', prob=TRUE, main='')
    lines(density(x$boot.est$par.boot[,par]))
  }else{
    hist(x$boot.est$boot.fit[,par], xlab=par, ylab='Density', prob=TRUE,main='')
    lines(density(x$boot.est$boot.fit[,par]))
  }	
}

plot.bmem<-function(x, par, ...){
  ## x: output from bmem function
  ## par: a parameter or a fit indice to use
  parnames<-colnames(x$boot.est$par.boot)
  fitnames<-colnames(x$boot.est$boot.fit)
  
  if (par %in% parnames){
    hist(x$boot.est$par.boot[,par], xlab=par, ylab='Density', prob=TRUE, main='')
    lines(density(x$boot.est$par.boot[,par]))
  }else{
    hist(x$boot.est$boot.fit[,par], xlab=par, ylab='Density', prob=TRUE,main='')
    lines(density(x$boot.est$boot.fit[,par]))
  }	
}

bmem.cov <- function(ram,obs.variables, debug=FALSE){
  ##if (moment){
  ##fixed.x<-'(i)'
  #3}else{
  fixed.x=NULL
  ##}
  parse.path <- function(path) {                                           
    path.1 <- gsub('-', '', gsub(' ','', path))
    direction <- if (regexpr('<>', path.1) > 0) 2 
    else if (regexpr('<', path.1) > 0) -1
    else if (regexpr('>', path.1) > 0) 1
    else stop(paste('ill-formed path:', path))
    path.1 <- strsplit(path.1, '[<>]')[[1]]
    list(first=path.1[1], second=path.1[length(path.1)], direction=direction)
  }
  if ((!is.matrix(ram)) | ncol(ram) != 3) stop ('ram argument must be a 3-column matrix')
  startvalues <- as.numeric(ram[,3])
  par.names <- ram[,2]
  n.paths <- length(par.names)
  heads <- from <- to <- rep(0, n.paths)
  for (p in 1:n.paths){
    path <- parse.path(ram[p,1])
    heads[p] <- abs(path$direction)
    to[p] <- path$second
    from[p] <- path$first
    if (path$direction == -1) {
      to[p] <- path$first
      from[p] <- path$second
    }
  }
  ram <- matrix(0, p, 5)
  all.vars <- unique(c(to, from))
  latent.vars <- setdiff(all.vars, obs.variables)
  not.used <- setdiff(obs.variables, all.vars)
  
  vars <- c(obs.variables, latent.vars)
  pars <- na.omit(unique(par.names))
  ram[,1] <- heads
  ram[,2] <- apply(outer(vars, to, '=='), 2, which)
  ram[,3] <- apply(outer(vars, from, '=='), 2, which)   
  par.nos <- apply(outer(pars, par.names, '=='), 2, which)
  if (length(par.nos) > 0)
    ram[,4] <- unlist(lapply(par.nos, function(x) if (length(x) == 0) 0 else x))
  ram[,5]<- startvalues
  colnames(ram) <- c('heads', 'to', 'from', 'parameter', 'start')
  if (!is.null(fixed.x)) fixed.x <- apply(outer(vars, fixed.x, '=='), 2, which)
  n <- length(obs.variables)
  m <- length(all.vars)
  t <- length(pars)
  if (debug) {
    cat('\n observed variables:\n') 
    print(paste(paste(1:n,':', sep=''), obs.variables, sep=''))
    cat('\n')
    if (m > n){ 
      cat('\n latent variables:\n')
      print(paste(paste((n+1):m,':', sep=''), latent.vars, sep=''))
      cat('\n')
    }
    cat('\n parameters:\n') 
    print(paste(paste(1:t,':', sep=''), pars, sep=''))
    cat('\n\n RAM:\n')
    print(ram)
  }
  n<-length(obs.variables)
  observed <- 1:n
  n.fix <- length(fixed.x)
  
  m <- max(ram[,c(2,3)])
  missing.variances <- setdiff(1:m, ram[,2][ram[,2] == ram[,3]])
  if (length(missing.variances) > 0) warning(paste( "\nThe model is almost surely misspecified; check also for missing covariances.\n"))
  t <- max(ram[,4])
  df <- n*(n + 1)/2 - t - n.fix*(n.fix + 1)/2
  if (df < 0) stop(paste("The model has negative degrees of freedom =", df))
  J <- matrix(0, n, m)
  correct <- matrix(2, m, m)
  diag(correct) <- 1
  J[cbind(1:n, observed)]<-1
  par.posn <-  sapply(1:t, function(i) which(ram[,4] == i)[1])
  colnames(ram)<-c("heads", "to", "from", "parameter", "start value")
  rownames(ram)<-rep("",nrow(ram))
  #if (length(param.names) > 0) rownames(ram)[par.posn]<-param.names
  fixed <- ram[,4] == 0
  sel.free <- ram[,4]
  sel.free[fixed] <- 1
  one.head <- ram[,1] == 1
  one.free <- which( (!fixed) & one.head )
  two.free <- which( (!fixed) & (!one.head) )
  arrows.1 <- ram[one.head, c(2,3), drop=FALSE]
  arrows.2 <- ram[!one.head, c(2,3), drop=FALSE]
  arrows.2t <- ram[!one.head, c(3,2), drop=FALSE]
  arrows.1.free <- ram[one.free,c(2,3), drop=FALSE]
  arrows.2.free <- ram[two.free,c(2,3), drop=FALSE]
  sel.free.1 <- sel.free[one.free]
  sel.free.2 <- sel.free[two.free]
  unique.free.1 <- unique(sel.free.1)
  unique.free.2 <- unique(sel.free.2)
  
  A <- P <- matrix(0, m, m)
  val <- ifelse (fixed, ram[,5], 0)
  A[arrows.1] <- val[one.head]
  P[arrows.2t] <- P[arrows.2] <- val[!one.head]
  I.Ainv <- solve(diag(m) - A)
  C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
  rownames(C)<-colnames(C)<-obs.variables
  C
}


summary.bmem <- function(object, boot.cl=TRUE, estimates=TRUE,...){
  par <- object$ci     ##list: (par), estimate, se.boot, 2.5%, 97.5%
  fit <- object$ci.fit ##list: (sta), estimate, se.boot, 2.5%, 97.5%
  infor <- object$infor ##list: samplesize, method, boot, sucboot, ci
  lci <- par[,"2.5%"]
  rci <- par[,"97.5%"]
  npar <- nrow(par)
  nfit <- nrow(fit)
  parname <- names(par[,1])
  staname <- names(fit[,1])
  lavpartable <- object$lavpartable
  est <- lavpartable$est   ##lavpartable:id,lhs,op,rhs,user,block,group,free,ustart,exo,label,plabel,start,est,se
  se  <- lavpartable$se
  
  cat("\n")
  t0.txt <- sprintf("Estimate method:                          ")
  t1.txt <- sprintf("%s",infor$method)
  cat(t0.txt,t1.txt, "\n", sep="")
  
  t0.txt <- sprintf("Sample size:                              ")
  t1.txt <- sprintf("%-10.d",infor$samplesize)
  cat(t0.txt,t1.txt, "\n", sep="")
  
  t0.txt <- sprintf("Number of request bootstrap draws:        ")
  t1.txt <- sprintf("%-10.d",infor$boot)
  cat(t0.txt,t1.txt, "\n", sep="")
  
  t0.txt <- sprintf("Number of successful bootstrap draws:     ")
  t1.txt <- sprintf("%-10.d",infor$sucboot)
  cat(t0.txt,t1.txt, "\n", sep="")
  
  t0.txt <- sprintf("Type of confidence interval:              ")
  t1.txt <- sprintf("%s",infor$ci)
  cat(t0.txt,t1.txt, "\n", sep="")
  
  cat("\n")
  
  
  if(boot.cl==FALSE){
    
    t0.txt <- sprintf("  %-40s", "Parameter")
    t1.txt <- sprintf("  %10s", "Estimate")
    cat(t0.txt, t1.txt, "\n", sep="")
    
    for (i in 1:npar) {
      t0.txt <- sprintf("  %-40s", parname[i])  
      t1.txt <- sprintf("  %10.3f", par[i,1])
      cat(t0.txt, t1.txt, "\n", sep="")
    }
    
    cat("\n")
    
    t0.txt <- sprintf("  %-40s", "Statistics")
    t1.txt <- sprintf("  %10s", "Estimate")
    cat(t0.txt, t1.txt, "\n", sep="")
    
    for (i in 1:nfit) {
      t0.txt <- sprintf("  %-40s", staname[i])  
      t1.txt <- sprintf("  %10.3f", fit[i,1])
      cat(t0.txt, t1.txt, "\n", sep="")
    }
  }
  
  if(boot.cl==TRUE){
    
    lpro <- colnames(par)[3]
    rpro <- colnames(par)[4]  ##confidence level 
    t0.txt <- sprintf("confidence level: %s ", rpro)
    cat(t0.txt, "\n", sep="")
    cat("\n")
    
    t0.txt <- sprintf("  %-20s", "Parameter")
    t1.txt <- sprintf("  %-20s", "Estimate")
    t2.txt <- sprintf("  %s", "boot.cl")
    cat(t0.txt, t1.txt, t2.txt, "\n", sep="")
    
    for (i in 1:npar) {
      t0.txt <- sprintf("  %-20s", parname[i])  
      t1.txt <- sprintf("  %-20.3f", par[i,1])
      t2.txt <- sprintf("  (%.3f, %.3f)", par[i,3],par[i,4]) ##stop here
      cat(t0.txt, t1.txt, t2.txt, "\n", sep="")
    }
  }
  
  cat("\n")
  t0.txt <- sprintf("Values of statistics:")
  cat(t0.txt, "\n", sep="")
  cat("\n")
  
  cat("                     Value      SE      ")
  cat(lpro)
  cat("     ")
  cat(rpro)
  cat("\n")
  for (i in 1:nfit) {
    t0.txt <- sprintf("  %-20s", staname[i])  
    t1.txt <- sprintf("%-8.3f", fit[i,1])
    t2.txt <- sprintf(" %-8.3f", fit[i,2])
    t3.txt <- sprintf(" %-8.3f", fit[i,3])
    t4.txt <- sprintf(" %-8.3f", fit[i,4])
    cat(t0.txt, t1.txt, t2.txt, t3.txt, t4.txt, "\n", sep="")
  }
  cat("\n")
  
  standardized <- 0
  
  if(estimates) {       
    t0.txt <- sprintf("Estimation of parameters:")
    cat(t0.txt, "\n", sep="")
    cat("\n")
    
    # local print function
    print.estimate <- function(name="ERROR", i=1, z.stat=TRUE) { ##name is the output of function makename
      
      # cut name if (still) too long
      name <- substr(name, 1, 13)
      
      if(!standardized) {

        if(is.na(se[i])) { ##
          txt <- sprintf("    %-13s %9.3f %8.3f\n", name, est[i], se[i])
        } else if(se[i] == 0) {
          txt <- sprintf("    %-13s %9.3f\n", name, est[i])
        } else if(est[i]/se[i] > 9999.999) {
          txt <- sprintf("    %-13s %9.3f %8.3f\n", name, est[i], se[i])
        } else if(!z.stat) {
          txt <- sprintf("    %-13s %9.3f %8.3f\n", name, est[i], se[i])
        } else {
          ##z <- est[i]/se[i]
          ##pval <- 2 * (1 - pnorm( abs(z) ))
          txt <- sprintf("    %-13s %9.3f %8.3f %8.3f %8.3f\n",
                         name, est[i], se[i], lci[i], rci[i])
        }
      } 
      cat(txt)
    } ##print.estimate end
    
    # estimates header
    if(!standardized) {
      cat("                     Estimate   SE      ")
      cat(lpro)
      cat("     ")
      cat(rpro)
      cat("\n")
    }
    
    ngroups <- length(unique(lavpartable$group))
    for(g in 1:ngroups) {
      
      makeNames <- function(NAMES, LABELS) {
        # labels?
        l.idx <- which(nchar(LABELS) > 0L)   ##select items with labels
        if(length(l.idx) > 0L) {
          LABELS <- abbreviate(LABELS, 4)    ##name the items
          LABELS[l.idx] <- paste(" (", LABELS[l.idx], ")", sep="") ##add a "(" ")"
          MAX.L <- max(nchar(LABELS))
          NAMES <- abbreviate(NAMES, minlength = (13 - MAX.L), strict = TRUE) ##name the items
          NAMES <- sprintf(paste("%-", (13 - MAX.L), "s%", MAX.L, "s", sep=""), NAMES, LABELS)
        } else {
          NAMES <- abbreviate(NAMES, minlength = 13, strict = TRUE)
        }
      }
      
      NAMES <- lavpartable$rhs
      
      # 1a. indicators ("=~") (we do show dummy indicators)
      mm.idx <- which( lavpartable$op == "=~" & 
                         ##!lavpartable$lhs %in% ov.names &
                         lavpartable$group == g)
      if(length(mm.idx)) {
        cat("Latent variables:\n")
        lhs.old <- ""
        NAMES[mm.idx] <- makeNames(  lavpartable$rhs[mm.idx],
                                     lavpartable$label[mm.idx] )
        for(i in mm.idx) {
          lhs <- lavpartable$lhs[i]
          if(lhs != lhs.old) cat("  ", lhs, " =~\n", sep="")
          print.estimate(name=NAMES[i], i)
          lhs.old <- lhs
        }
        cat("\n")
      }
      
      # 1b. formative/composites ("<~")
      fm.idx <- which( lavpartable$op == "<~" & lavpartable$group == g)
      if(length(fm.idx)) {
        cat("Composites:\n")
        lhs.old <- ""
        NAMES[fm.idx] <- makeNames(lavpartable$rhs[fm.idx], lavpartable$label[fm.idx])
        for(i in fm.idx) {
          lhs <- lavpartable$lhs[i]
          if(lhs != lhs.old) cat("  ", lhs, " <~\n", sep="")
          print.estimate(name=NAMES[i], i)
          lhs.old <- lhs
        }
        cat("\n")
      }
      
      # 2. regressions
      eqs.idx <- which(lavpartable$op == "~" & lavpartable$group == g)
      if(length(eqs.idx) > 0) {
        cat("Regressions:\n")
        lhs.old <- ""
        NAMES[eqs.idx] <- makeNames(  lavpartable$rhs[eqs.idx],
                                      lavpartable$label[eqs.idx])
        for(i in eqs.idx) {
          lhs <- lavpartable$lhs[i]
          if(lhs != lhs.old) cat("  ", lhs, " ~\n", sep="")
          print.estimate(name=NAMES[i], i)
          lhs.old <- lhs
        }
        cat("\n")
      }
      
      # 3. covariances
      cov.idx <- which(lavpartable$op == "~~" & 
                         ##!lavpartable$exo &
                         lavpartable$lhs != lavpartable$rhs &
                         lavpartable$group == g)
      if(length(cov.idx) > 0) {
        cat("Covariances:\n")
        lhs.old <- ""
        NAMES[cov.idx] <- makeNames(  lavpartable$rhs[cov.idx],
                                      lavpartable$label[cov.idx])
        for(i in cov.idx) {
          lhs <- lavpartable$lhs[i]
          if(lhs != lhs.old) cat("  ", lhs, " ~~\n", sep="")
          print.estimate(name=NAMES[i], i)
          lhs.old <- lhs
        }
        cat("\n")
      }
      
      # 4. intercepts/means
      int.idx <- which(lavpartable$op == "~1" & 
                         ##!lavpartable$exo &
                         lavpartable$group == g)
      if(length(int.idx) > 0) {
        cat("Intercepts:\n")
        NAMES[int.idx] <- makeNames(  lavpartable$lhs[int.idx],
                                      lavpartable$label[int.idx])
        for(i in int.idx) {
          print.estimate(name=NAMES[i], i)
        }
        cat("\n")
      }
      
      # 5. (residual) variances
      var.idx <- which(lavpartable$op == "~~" &
                         ##!lavpartable$exo &
                         lavpartable$lhs == lavpartable$rhs &
                         lavpartable$group == g)
      if(length(var.idx) > 0) {
        cat("Variances:\n")
        NAMES[var.idx] <- makeNames(  lavpartable$rhs[var.idx],
                                      lavpartable$label[var.idx])
        for(i in var.idx) {
          ##if(lavoptions$mimic == "lavaan") {
          ##print.estimate(name=NAMES[i], i, z.stat=TRUE)
          ##} else {
          print.estimate(name=NAMES[i], i, z.stat=TRUE)
        }
      }
      cat("\n")
    } # "for" of ngroups end
    
    
    # 6. variable definitions
    def.idx <- which(lavpartable$op == ":=")
    if(length(def.idx) > 0) {
      if(ngroups > 1) cat("\n")
      cat("Defined parameters:\n")
      ##NAMES[def.idx] <- makeNames(  lavpartable$lhs[def.idx], "")
      NAMES[def.idx] <- makeNames(  lavpartable$rhs[def.idx], lavpartable$label[def.idx])
      for(i in def.idx) {
        print.estimate(name=NAMES[i], i)
      }
      cat("\n")
    }
    
    # 7. constraints
    cin.idx <- which((lavpartable$op == "<" | 
                        lavpartable$op == ">"))
    ceq.idx <- which(lavpartable$op == "==")
    if(length(cin.idx) > 0L || length(ceq.idx) > 0L) {
      # set small negative values to zero, to avoid printing " -0.000"
      slack <- ifelse(abs(est) < 1e-5, 0, est)
      
      if(ngroups > 1 && length(def.idx) == 0L) cat("\n")
      cat("Constraints:                               Slack (>=0)\n")
      for(i in c(cin.idx,ceq.idx)) {
        lhs <- lavpartable$lhs[i]
        op <- lavpartable$op[i]
        rhs <- lavpartable$rhs[i]
        if(rhs == "0" && op == ">") {
          con.string <- paste(lhs, " - 0", sep="")
        } else if(rhs == "0" && op == "<") {
          con.string <- paste(rhs, " - (", lhs, ")", sep="")
        } else if(rhs != "0" && op == ">") {
          con.string <- paste(lhs, " - (", rhs, ")", sep="")
        } else if(rhs != "0" && op == "<") {
          con.string <- paste(rhs, " - (", lhs, ")", sep="")
        } else if(rhs == "0" && op == "==") {
          con.string <- paste(lhs, " - 0", sep="")
        } else if(rhs != "0" && op == "==") {
          con.string <- paste(lhs, " - (", rhs, ")", sep="")
        }
        con.string <- abbreviate(con.string, 41, strict = TRUE)
        txt <- sprintf("    %-41s %8.3f\n", 
                       con.string, slack[i])
        cat(txt)
      }   
      cat("\n")
    }
    
  } # "if(estimates)" end
  
}