====== Some results ======
<code>structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.059, 0.072, 0.058, 
0.046, 0.057, 0.065, 0.069, 0.043, 0.055, 0.054, 0.072, 0.059, 
0.058, 0.05, 0.062, 0.059, 0.046, 0.073, 0.042, 0.051, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0.972, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.814, 
0.974, 1, 1, 1, 1, 1, 1, 1, 1, 0.146, 0.206, 0.288, 0.384, 0.452, 
0.528, 0.593, 0.651, 0.69, 0.746, 0.291, 0.505, 0.665, 0.794, 
0.897, 0.928, 0.952, 0.979, 0.987, 0.991, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 0.442, 0.714, 0.896, 0.962, 0.976, 0.991, 0.997, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0.105, 0.188, 0.239, 0.33, 0.395, 0.474, 
0.499, 0.544, 0.598, 0.65, 0.135, 0.2, 0.236, 0.318, 0.366, 0.432, 
0.484, 0.52, 0.607, 0.622), .Dim = c(10L, 21L))</code>

====== R code ======
<code>

## Monte Calor based power analysis for indirect effect using bootstrap
## by Zhiyong Zhang

## bootstrap confidence intervals
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))
	}

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
	}

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
	}

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, ci.bc1(par.boot[,i], par0[i], cl))
		}
	cbind(estimate, se.boot, ci)
	}

ci.sig<-function(res){
	sig<-function(a){
		ifelse(a[1]>0 | a[2]<0, 1, 0)
	}
	apply(res[, 3:4], 1, sig)
}

popPar<-function(object){
	par.value<-object@ParTable$ustart
	for (i in 1:length(par.value)){
		if (is.na(par.value[i])){
			if (object@ParTable$op[i]=="~"){
				par.value[i]<-0
			}else{
				par.value[i]<-1
			}
		}
	}
	names(par.value)<-object@ParTable$label
	## for indirect effec defined here
	for (i in 1:length(par.value)){
		if (object@ParTable$op[i]==":="){
			temp<-object@ParTable$rhs[i]
			temp<-gsub(' +','',temp)
			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=object@ParTable$rhs[i])
			par.list<-as.list(par.value[par])
			par.value[i]<-eval(ind.exp, par.list)
		}
	}
	par.value
}

## power function for analysis
power.basic<-function(model, indirect=NULL, nobs, nrep=1000, alpha=.95, skewness=NULL, kurtosis=NULL, se="default", estimator="default", ...){
	if (missing(model)) stop("A model is needed.")
	## repeated the following for nrep times
	all.sig<-NULL
	all.par<-NULL
	all.se<-NULL
	all.cvg<-NULL
	
	model.indirect<-paste(model, "\n", indirect, "\n")
	## Initial analysis for some model information
	newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)			
	temp.res<-sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
	par.value<-popPar(temp.res)
	
	for (i in 1:nrep){
		## Step 1: generate data
		newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)		
		## Step 2: fit the model 		
		temp.res<-sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
		
		## Step 3: Check significance
		idx <- 1:length( temp.res@ParTable$lhs )
		temp.est<-temp.res@Fit@est[idx]
		temp.se<-temp.res@Fit@se[idx]
		temp.sig<-temp.est/temp.se
		crit<-qnorm(1-(1-alpha)/2)
		temp.sig.ind<-abs(temp.sig)>crit
		temp.sig.ind[!is.finite(temp.sig)]<-NA
		
		## Step 4: Check the coverage
		ci.u<-temp.est+crit*temp.se
		ci.l<-temp.est-crit*temp.se
		temp.cvg<- (ci.u>par.value[idx]) & (ci.l<par.value[idx])
		
		## Step 5: Combine results
		all.sig<-rbind(all.sig, temp.sig.ind)
		all.par<-rbind(all.par, temp.est)
		all.se<-rbind(all.se, temp.se)
		all.cvg<-rbind(all.cvg, temp.cvg)
	}
	cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
	colnames(all.sig)<-cnames
	power<-apply(all.sig, 2, mean, na.rm=TRUE)
	cvg<-apply(all.cvg, 2, mean, na.rm=TRUE)
	info<-list(nobs=nobs, nrep=nrep, alpha=alpha, method="Normal", bootstrap=NULL)
	print(power)
	object<-list(power=power, coverage=cvg, pop.value=par.value, results=list(estimates=all.par, se=all.se), info=info, out=temp.res)
	class(object)<-'power'
	return(object)
}

power.boot<-function(model, indirect=NULL, nobs, nrep=1000, nboot=1000, alpha=.95, skewness=NULL, kurtosis=NULL, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), ...){		
	if (missing(model)) stop("A model is needed.")	
	all.sig<-NULL
	all.par<-NULL
	all.se<-NULL
	all.cvg<-NULL
	ci.all<-list()
	boot.all<-list()
	## internal function
	coef.new<-function(x,...){
		coef(x, ...)
	}
	model.indirect<-paste(model, "\n", indirect, "\n")
	## Initial analysis for some model information
	newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)			
	temp.res<-sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
	par.value<-popPar(temp.res)
	idx <- 1:length( temp.res@ParTable$lhs )
	
	## repeated the following for nrep times
	for (i in 1:nrep){
		## Step 1: generate data
		newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)	
		
		## Step 2: fit the model 		
		temp.res<-sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
		
		## Step 3: Conduct bootstrap analysis
		orig.res<-coef.new(temp.res)
		boot.res<-bootstrapLavaan(temp.res, FUN=coef.new, R=nboot, parallel=parallel, ncpus=ncore, ...)
		ci.res<-ci.bc(boot.res, orig.res, cl=alpha)
		ci.all[[i]]<-ci.res
		boot.all[[i]]<-boot.res
		## Step 3: Check significance		
		temp.est<-temp.res@Fit@est[idx]
		temp.se<-ci.res[,2]
		temp.sig<-ci.sig(ci.res)
		
		## Step 4: Check the coverage
		temp.cvg<- (ci.res[,4]>par.value[idx]) & (ci.res[,3]<par.value[idx])
		
		## Step 5: Combine results
		all.sig<-rbind(all.sig, temp.sig)
		all.par<-rbind(all.par, temp.est)
		all.se<-rbind(all.se, temp.se)
		all.cvg<-rbind(all.cvg, temp.cvg)
	}
	cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
	colnames(all.sig)<-cnames
	power<-apply(all.sig, 2, mean, na.rm=TRUE)
	cvg<-apply(all.cvg, 2, mean, na.rm=TRUE)
	info<-list(nobs=nobs, nrep=nrep, alpha=alpha, method="Normal", bootstrap=nboot)
	print(power)
	object<-list(power=power, coverage=cvg, pop.value=par.value, results=list(estimates=all.par, se=all.se, ci=ci.all, boot=boot.all), info=info, out=temp.res)
	class(object)<-'power'
	return(object)
}



summary.power<-function(object) {

    # Show some basic about the simulation methods

    # main part: parameter estimates
    cat("Basic information:\n\n")
	if(object$out@Options$test %in% c("satorra.bentler", "yuan.bentler",
                                  "mean.var.adjusted",
                                  "scaled.shifted") &&
       length(object$out@Fit@test) > 1L) {
        scaled <- TRUE
        if(object$out@Options$test == "scaled.shifted")
            shifted <- TRUE
        else
            shifted <- FALSE
    } else {
        scaled <- FALSE
        shifted <- FALSE
    }
	
	t0.txt <- sprintf("  %-40s", "Esimation method")
    t1.txt <- sprintf("  %10s", object$out@Options$estimator)
    t2.txt <- ifelse(scaled, 
              sprintf("  %10s", "Robust"), "")
    cat(t0.txt, t1.txt, t2.txt, "\n", sep="")
	
	t0.txt <- sprintf("  %-40s", "Standard error")
    t1.txt <- sprintf("  %10s", object$out@Options$se)
	cat(t0.txt, t1.txt, "\n", sep="")
	
	if (!is.null(object$info$bootstrap)){
		t0.txt <- sprintf("  %-40s", "Number of requested bootstrap")
		t1.txt <- sprintf("  %10i", object$info$bootstrap)
		cat(t0.txt, t1.txt, "\n", sep="")
	}

    t0.txt <- sprintf("  %-40s", "Number of requested replications")
    t1.txt <- sprintf("  %10i", object$info$nrep)
    cat(t0.txt, t1.txt, "\n", sep="")
    t0.txt <- sprintf("  %-40s", "Number of successful replications")
    t1.txt <- sprintf("  %10i", nrow(object$results$estimates))
    cat(t0.txt, t1.txt, "\n", sep="")
		
    cat("\n")

    # local print function
    print.estimate <- function(name="ERROR", i=1, z.stat=TRUE) {       
        # cut name if (still) too long
        name <- strtrim(name, width=15L)
        txt <- sprintf("    %-15s %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
                               name, pop.value[i], est[i], mse[i], sd.est[i], power[i], cvg[i])
        cat(txt)
    }

    est <- apply(object$results$estimates, 2, mean, na.rm=TRUE)
    sd.est  <- apply(object$results$estimates, 2, sd, na.rm=TRUE)
	mse <- apply(object$results$se, 2, mean, na.rm=TRUE)
	power <- object$power
	pop.value<-object$pop.value
	cvg<-object$coverage
	
	object<-object$out
	

    for(g in 1:object@Data@ngroups) {
        ov.names <- lavaan:::vnames(object@ParTable, "ov", group=g)
        lv.names <- lavaan:::vnames(object@ParTable, "lv", group=g)

        # group header
        if(object@Data@ngroups > 1) {
            if(g > 1) cat("\n\n")
            cat("Group ", g, 
                " [", object@Data@group.label[[g]], "]:\n\n", sep="")
        }

        # estimates header
		#txt <- sprintf("%-13s %12s %12s %8s %8s %8s\n", "", "True", "Estimate", "MSE", "SD", "Power")
		#cat(txt)
        cat("                       True  Estimate      MSE      SD     Power Coverage\n")
		
        makeNames <- function(NAMES, LABELS) {
            multiB <- FALSE
            if(any(nchar(NAMES) != nchar(NAMES, "bytes")))
                multiB <- TRUE
            if(any(nchar(LABELS) != nchar(LABELS, "bytes")))
                multiB <- TRUE
            # labels?
            l.idx <- which(nchar(LABELS) > 0L)
            if(length(l.idx) > 0L) {
                if(!multiB) {
                    LABELS <- abbreviate(LABELS, 4)
                    LABELS[l.idx] <- paste(" (", LABELS[l.idx], ")", sep="")
                    MAX.L <- max(nchar(LABELS))
                    NAMES <- abbreviate(NAMES, minlength = (13 - MAX.L), 
                                        strict = TRUE)
                } else {
                    # do not abbreviate anything (eg in multi-byte locales)
                    MAX.L <- 4L
                }
                NAMES <- sprintf(paste("%-", (13 - MAX.L), "s%", MAX.L, "s",
                                       sep=""), NAMES, LABELS)
            } else {
                if(!multiB) {
                    NAMES <- abbreviate(NAMES, minlength = 13, strict = TRUE)
                } else {
                    NAMES <- sprintf(paste("%-", 13, "s", sep=""), NAMES)
                }
            }

            NAMES
        }

        NAMES <- object@ParTable$rhs

        # 1a. indicators ("=~") (we do show dummy indicators)
        mm.idx <- which( object@ParTable$op == "=~" & 
                        !object@ParTable$lhs %in% ov.names &
                         object@ParTable$group == g)
        if(length(mm.idx)) {
            cat("Latent variables:\n")
            lhs.old <- ""
            NAMES[mm.idx] <- makeNames(  object@ParTable$rhs[mm.idx],
                                       object@ParTable$label[mm.idx])
            for(i in mm.idx) {
                lhs <- object@ParTable$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( object@ParTable$op == "<~" &
                         object@ParTable$group == g)
        if(length(fm.idx)) {
            cat("Composites:\n")
            lhs.old <- ""
            NAMES[fm.idx] <- makeNames(  object@ParTable$rhs[fm.idx],
                                       object@ParTable$label[fm.idx])
            for(i in fm.idx) {
                lhs <- object@ParTable$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(object@ParTable$op == "~" & object@ParTable$group == g)
        if(length(eqs.idx) > 0) {
            cat("Regressions:\n")
            lhs.old <- ""
            NAMES[eqs.idx] <- makeNames(  object@ParTable$rhs[eqs.idx],
                                        object@ParTable$label[eqs.idx])
            for(i in eqs.idx) {
                lhs <- object@ParTable$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(object@ParTable$op == "~~" & 
                         !object@ParTable$exo &
                         object@ParTable$lhs != object@ParTable$rhs &
                         object@ParTable$group == g)
        if(length(cov.idx) > 0) {
            cat("Covariances:\n")
            lhs.old <- ""
            NAMES[cov.idx] <- makeNames(  object@ParTable$rhs[cov.idx],
                                        object@ParTable$label[cov.idx])
            for(i in cov.idx) {
                lhs <- object@ParTable$lhs[i]
                if(lhs != lhs.old) cat("  ", lhs, " ~~\n", sep="")
                print.estimate(name=NAMES[i], i)
                lhs.old <- lhs
            }
            cat("\n")
        }

        # 4. intercepts/means
        ord.names <- lavaan:::vnames(object@ParTable, type="ov.ord", group=g)
        int.idx <- which(object@ParTable$op == "~1" & 
                         !object@ParTable$lhs %in% ord.names &
                         !object@ParTable$exo &
                         object@ParTable$group == g)
        if(length(int.idx) > 0) {
            cat("Intercepts:\n")
            NAMES[int.idx] <- makeNames(  object@ParTable$lhs[int.idx],
                                        object@ParTable$label[int.idx])
            for(i in int.idx) {
                print.estimate(name=NAMES[i], i)
            }
            cat("\n")
        }

        # 4b thresholds
        th.idx <- which(object@ParTable$op == "|" &
                        object@ParTable$group == g)
        if(length(th.idx) > 0) {
            cat("Thresholds:\n")
            NAMES[th.idx] <- makeNames(  paste(object@ParTable$lhs[th.idx],
                                               "|",
                                               object@ParTable$rhs[th.idx],
                                               sep=""),
                                         object@ParTable$label[th.idx])
            for(i in th.idx) {
                print.estimate(name=NAMES[i], i)
            }
            cat("\n")
        }

        # 5. (residual) variances
        var.idx <- which(object@ParTable$op == "~~" &
                         !object@ParTable$exo &
                         object@ParTable$lhs == object@ParTable$rhs &
                         object@ParTable$group == g)
        if(length(var.idx) > 0) {
            cat("Variances:\n")
            NAMES[var.idx] <- makeNames(  object@ParTable$rhs[var.idx],
                                        object@ParTable$label[var.idx])
            for(i in var.idx) {
                if(object@Options$mimic == "lavaan") {
                    print.estimate(name=NAMES[i], i, z.stat=FALSE)
                } else {
                    print.estimate(name=NAMES[i], i, z.stat=TRUE)
                }
            }
            cat("\n")
        }

        # 6. latent response scales
        delta.idx <- which(object@ParTable$op == "~*~" &
                         object@ParTable$group == g)
        if(length(delta.idx) > 0) {
            cat("Scales y*:\n")
            NAMES[delta.idx] <- makeNames(  object@ParTable$rhs[delta.idx],
                                            object@ParTable$label[delta.idx])
            for(i in delta.idx) {
                print.estimate(name=NAMES[i], i, z.stat=TRUE)
            }
            cat("\n")
        }

    } # ngroups

    # 6. variable definitions
    def.idx <- which(object@ParTable$op == ":=")
    if(length(def.idx) > 0) {
        if(object@Data@ngroups > 1) cat("\n")
        cat("Indirect effects:\n")
        NAMES[def.idx] <- makeNames(  object@ParTable$lhs[def.idx], "")
        for(i in def.idx) {
            print.estimate(name=NAMES[i], i)
        }
        cat("\n")
    }

}

</code>