R/power.R

In bmem: Mediation Analysis with Missing Data Using Bootstrap

## Monte Calor based power analysis for indirect effect using bootstrap
## by Zhiyong Zhang
summary.power<-function(object, ...) {

    # lavaan slots
    lavoptions <- lavInspect(object$out, "options")
    lavtest    <- lavInspect(object$out, "test")

    # Show some basic about the simulation methods

    # main part: parameter estimates
    cat("Basic information:\n\n")
	if(lavoptions$test %in% c("satorra.bentler", "yuan.bentler",
                              "mean.var.adjusted",
                              "scaled.shifted") &&
       length(lavtest) > 1L) {
        scaled <- TRUE
        if(lavoptions$test == "scaled.shifted")
            shifted <- TRUE
        else
            shifted <- FALSE
    } else {
        scaled <- FALSE
        shifted <- FALSE
    }
	
	t0.txt <- sprintf("  %-40s", "Esimation method")
    t1.txt <- sprintf("  %10s", lavoptions$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", lavoptions$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 %8.3f\n",
                               name, pop.value[i], est[i], mse[i], sd.est[i], power[i], power.se[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
	power.se<-object$power.se
	pop.value<-object$pop.value
	cvg<-object$coverage
	
	object<-object$out
	ngroups <- lavInspect(object, "ngroups")
    group.label <- lavInspect(object, "group.label")
    lavpartable <- parTable(object)

    for(g in 1:ngroups) {
        ov.names <- lavaanNames(object, group = g) 
        lv.names <- lavaanNames(object, type = "lv", group = g) 

        # group header
        if(ngroups > 1) {
            if(g > 1) cat("\n\n")
            cat("Group ", g, 
                " [", 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 Power.se 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 <- 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
        ord.names <- lavaanNames(object, type="ov.ord") 
        int.idx <- which(lavpartable$op == "~1" & 
                         !lavpartable$lhs %in% ord.names &
                         !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")
        }

        # 4b thresholds
        th.idx <- which(lavpartable$op == "|" &
                        lavpartable$group == g)
        if(length(th.idx) > 0) {
            cat("Thresholds:\n")
            NAMES[th.idx] <- makeNames(  paste(lavpartable$lhs[th.idx],
                                               "|",
                                               lavpartable$rhs[th.idx],
                                               sep=""),
                                               lavpartable$label[th.idx])
            for(i in th.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=FALSE)
                } else {
                    print.estimate(name=NAMES[i], i, z.stat=TRUE)
                }
            }
            cat("\n")
        }

        # 6. latent response scales
        delta.idx <- which(lavpartable$op == "~*~" &
                           lavpartable$group == g)
        if(length(delta.idx) > 0) {
            cat("Scales y*:\n")
            NAMES[delta.idx] <- makeNames(  lavpartable$rhs[delta.idx],
                                            lavpartable$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(lavpartable$op == ":=")
    if(length(def.idx) > 0) {
        if(ngroups > 1) cat("\n")
        cat("Indirect/Mediation effects:\n")
        NAMES[def.idx] <- makeNames(  lavpartable$lhs[def.idx], "")
        for(i in def.idx) {
            print.estimate(name=NAMES[i], i)
        }
        cat("\n")
    }

}

## bootstrap confidence intervals
popPar<-function(object){

    lavpartable <- parTable(object)

	par.value<-lavpartable$ustart
	for (i in 1:length(par.value)){
		if (is.na(par.value[i])){
			if (lavpartable$op[i]=="~~"){
				if (lavpartable$rhs[i]==lavpartable$lhs[i]){
					par.value[i]<-1
				}else{
					par.value[i]<-0
				}
			}else{
				par.value[i]<-0
			}
		}
	}
	names(par.value)<-lavpartable$label
	## for indirect effec defined here
	for (i in 1:length(par.value)){
		if (lavpartable$op[i]==":="){
			temp<-lavpartable$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=lavpartable$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=100, nrep=1000, alpha=.95, skewness=NULL, kurtosis=NULL, ovnames=NULL, se="default", estimator="default", parallel="no", ncore=1, ...){
	if (missing(model)) stop("A model is needed.")
	
	model.indirect<-paste(model, "\n", indirect, "\n")
	ngroups <- length(nobs)
	## Initial analysis for some model information
	
	## checking for error
	error <- 0
	while (error == 0){
		newdata<-try(simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...))
		if (!inherits(newdata,  "try-error")) error <- 1
	}

	if (ngroups > 1){
		error <- 0
		while (error == 0){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, group='group', ...))
			if (!inherits(temp.res,  "try-error")) error <- 1
		}
	}else{
		error <- 0
		while (error == 0){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, ...))
			if (!inherits(temp.res,  "try-error")) error <- 1
		}
	}
	par.value<-popPar(temp.res)

    temp.res.lavpartable <- parTable(temp.res)
	idx <- 1:length( temp.res.lavpartable$lhs )
	#cnames<-paste(temp.res.lavpartable$lhs[idx], temp.res.lavpartable$op[idx], temp.res.lavpartable$rhs[idx])
	ov<-lavaanNames(temp.res,type="ov")
	## dealing with skewness and kurtosis
	if (length(skewness) > 1){
		if (length(skewness)!=length(ovnames)) stop("The number of skewness is not equal to the number observed variables")	
		index<-match(ov, ovnames)
		skewness<-skewness[index]
	}
	
	if (length(kurtosis) > 1){
		if (length(kurtosis)!=length(ovnames)) stop("The number of kurtosis is not equal to the number observed variables")	
		index<-match(ov, ovnames)
		kurtosis<-kurtosis[index]
	}	
	
	cnames<-paste(temp.res.lavpartable$lhs[idx], temp.res.lavpartable$op[idx], temp.res.lavpartable$rhs[idx])
	out<-temp.res
	runonce<-function(i){
		## Step 1: generate data
		error <- 0
		while (error == 0){
			newdata<-try(simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)	)
			if (!inherits(newdata,  "try-error")) error <- 1
		}
	
		## Step 2: fit the model 		
		if (ngroups > 1){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, group='group',  ...))
		}else{
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator,  ...))
		}
		
		## Step 3: Check significance
		
		if (!inherits(temp.res,  "try-error")){
            temp.res.lavpartable <- parTable(temp.res)
			idx <- 1:length( temp.res.lavpartable$lhs )
			temp.est<-temp.res.lavpartable$est[idx]
			temp.se<-temp.res.lavpartable$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])
		}else{
			nna<-length(idx)
			temp.est<-rep(NA, nna)
			temp.sig.ind<-rep(NA, nna)
			temp.se<-rep(NA, nna)
			temp.cvg<-rep(NA, nna)
		}		
		return(list(temp.sig.ind=temp.sig.ind, temp.est=temp.est, temp.se=temp.se, temp.cvg=temp.cvg))
	}
	
	## run parallel or not
	# this is from the boot function in package boot
	old_options <- options(); options(warn = -1)
    RR <- nrep
    res <- if (ncore > 1L && parallel != "no") {
		sfInit( parallel=TRUE, cpus=ncore )
		sfLibrary("lavaan", character.only = TRUE)
        sfLapply(seq_len(RR), runonce)
    } else lapply(seq_len(RR), runonce)
	
	all.sig<-do.call(rbind, lapply(res, "[[", 'temp.sig.ind'))

	all.par<-do.call(rbind, lapply(res, "[[", 'temp.est'))
	all.se<-do.call(rbind, lapply(res, "[[", 'temp.se'))
	all.cvg<-do.call(rbind, lapply(res, "[[", 'temp.cvg'))
		
	options(old_options)	
	
	colnames(all.sig)<-cnames
	power<-apply(all.sig, 2, mean, na.rm=TRUE)
	power.se<-sqrt(power*(1-power)/apply(!is.na(all.sig), 2, sum))
	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, power.se=power.se, coverage=cvg, pop.value=par.value, results=list(estimates=all.par, se=all.se, all=res), info=info, out=out, data=newdata)
	class(object)<-'power'
	return(object)
}


power.boot<-function(model, indirect=NULL, nobs=100, nrep=1000, nboot=1000, alpha=.95, skewness=NULL, kurtosis=NULL, ovnames=NULL, ci='default', boot.type='default', se="default", estimator="default", parallel="no", ncore=1, ...){		
	if (missing(model)) stop("A model is needed.")	
	
	## internal function
	coef.new<-function(x,...){
		coef(x, type='user', ...)
	}
	model.indirect<-paste(model, "\n", indirect, "\n")
	ngroups <- length(nobs)
		
	## Initial analysis for some model information
	## checking for error
	error <- 0
	while (error == 0){
		newdata<-try(simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...))
		if (!inherits(newdata,  "try-error")) error <- 1
	}

	if (ngroups > 1){
		error <- 0
		while (error == 0){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, group='group', ...))
			if (!inherits(temp.res,  "try-error")) error <- 1
		}
	}else{
		error <- 0
		while (error == 0){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, ...))
			if (!inherits(temp.res,  "try-error")) error <- 1
		}
	}
    temp.res.lavpartable <- parTable(temp.res)
	par.value<-popPar(temp.res)
	idx <- 1:length( temp.res.lavpartable$lhs )
	#cnames<-paste(temp.res.lavpartable$lhs[idx], temp.res.lavpartable$op[idx], temp.res.lavpartable$rhs[idx])
	ov<-lavaanNames(temp.res,type="ov")
	ptype<-(temp.res.lavpartable$op == "~~") & (temp.res.lavpartable$lhs == temp.res.lavpartable$rhs)
	
	## dealing with skewness and kurtosis
	if (length(skewness) > 1){
		if (length(skewness)!=length(ovnames)) stop("The number of skewness is not equal to the number observed variables")	
		index<-match(ov, ovnames)
		skewness<-skewness[index]
	}
	
	if (length(kurtosis) > 1){
		if (length(kurtosis)!=length(ovnames)) stop("The number of kurtosis is not equal to the number observed variables")	
		index<-match(ov, ovnames)
		kurtosis<-kurtosis[index]
	}
		
	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.perc<-function(par.boot, par0, cl=.95){
		se.boot<-apply(par.boot, 2, sd, na.rm=TRUE)
		estimate<-par0
		alpha<-(1-cl)/2
		alpha<-c(alpha, 1-alpha)
		ci<-apply(par.boot, 2, quantile, prob=alpha) 
		cbind(estimate, se.boot, t(ci))
	}
	
	## repeated the following for nrep times
	runonce<-function(i){
		# Step 1: generate data
		error <- 0
		while (error == 0){
			newdata<-try(simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)	)
			if (!inherits(newdata,  "try-error")) error <- 1
		}
		
		## Step 2: fit the model 		
		#temp.res<-sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
		if (ngroups > 1){
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator, group='group',  ...))
		}else{
			temp.res<-try(lavaan::cfa(model.indirect, data=newdata, se=se, estimator=estimator,  ...))
		}
		
		
		## Step 3: Conduct bootstrap analysis
		if (!inherits(temp.res,  "try-error")){
			orig.res<-coef.new(temp.res)
			if (boot.type=="default"){
				boot.res<-bootstrapLavaan(temp.res, FUN=coef.new, R=nboot, parallel="no",  ...)
			}else{
				boot.res<-bootstrapLavaan(temp.res, FUN=coef.new, R=nboot, type='bollen.stine', parallel="no",  ...)
			}
			if (ci=='default'){
				ci.res<-ci.perc(boot.res, orig.res, cl=alpha)
			}else{
				ci.res<-ci.bc(boot.res, orig.res, cl=alpha)
			}
			## Step 4: Check the coverage		
		 
			temp.cvg<- (ci.res[,4]>=par.value[idx]) & (ci.res[,3]<par.value[idx])
		
			## Step 5: check significance
			temp.sig<-NULL
			temp.est<-coef.new(temp.res) 
			temp.se<-ci.res[,2]
			for (jj in 1:length(ptype)){
				if (ptype[jj]){
					crit<-qnorm(1-(1-alpha)/2)
					ci.temp<-ci.res[jj, 1] + c(-1,1)*ci.res[jj, 2]*crit
					temp.sig<-c(temp.sig, (ci.temp[1]>0 | ci.temp[2]<0))
				}else{
					temp.sig<-c(temp.sig, (ci.res[jj, 3]>0 | ci.res[jj, 4]<0))
				}
			}
		}else{
			nna<-length(idx)
			temp.est<-rep(NA, nna)
			temp.sig<-rep(NA, nna)
			temp.se<-rep(NA, nna)
			temp.cvg<-rep(NA, nna)
			boot.res<-NA
			ci.res<-NA
		}

		return(list(temp.sig=temp.sig, temp.est=temp.est, temp.se=temp.se, temp.cvg=temp.cvg, boot.res=boot.res, ci.res=ci.res))
	}
	
	## run parallel or not
	old_options <- options(); options(warn = -1)
	
    RR <- nrep
    res <- if (ncore > 1L && parallel != "no") {
		sfInit( parallel=TRUE, cpus=ncore )
		sfLibrary("lavaan", character.only = TRUE)
        sfLapply(seq_len(RR), runonce)
    } else lapply(seq_len(RR), runonce)
	
	all.sig<-do.call(rbind, lapply(res, "[[", 'temp.sig'))

	all.par<-do.call(rbind, lapply(res, "[[", 'temp.est'))
	all.se<-do.call(rbind, lapply(res, "[[", 'temp.se'))
	all.cvg<-do.call(rbind, lapply(res, "[[", 'temp.cvg'))
		
	options(old_options)	
	
	#colnames(all.sig)<-cnames
	power<-apply(all.sig, 2, mean, na.rm=TRUE)
	power.se<-sqrt(power*(1-power)/apply(!is.na(all.sig), 2, sum))
	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, power.se=power.se, coverage=cvg, pop.value=par.value, results=list(estimates=all.par, se=all.se, all=res), info=info, out=temp.res, data=newdata)
	class(object)<-'power'
	return(object)
}


power.curve<-function(model, indirect=NULL, nobs=100, type='basic', nrep=1000, nboot=1000, alpha=.95, skewness=NULL, kurtosis=NULL, ovnames=NULL,  ci='default', boot.type='default', se="default", estimator="default", parallel="no", ncore=1, interactive=TRUE, ...){		
	if (missing(model)) stop("A model is needed.")	
	
	allpower<-NULL
	
	## check whether nobs is a vector or not
	if (is.vector(nobs)){	
		for (N in nobs){
			if (type == 'basic'){
				## sobel test based analysis
				indpower <- power.basic(model=model, indirect=indirect, nobs=N, nrep=nrep, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore,  ...)			
			}else{
				indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, ci=ci, boot.type=boot.type, estimator=estimator, parallel=parallel, ncore=ncore,  ...)	
			}
			allpower<-rbind(allpower, indpower$power)
		}
	}else{
		for (k in 1:nrow(nobs)){
			N <- nobs[k]
			if (type == 'basic'){
				## sobel test based analysis
				indpower <- power.basic(model=model, indirect=indirect, nobs=N, nrep=nrep, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore,  ...)			
			}else{
				indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, ci=ci, boot.type=boot.type, estimator=estimator, parallel=parallel, ncore=ncore,  ...)	
			}
			allpower<-rbind(allpower, indpower$power)
		}
	}

	if (interactive & .Platform$OS.type=="windows"){
		dev.new(record=TRUE)
		op <- par(ask=TRUE)
		on.exit(par(op))
	}
	partable <- parTable(indpower$out)
	pnames <- paste0(partable$lhs, partable$op,  partable$rhs, " ", partable$label)
	for (j in 1:ncol(allpower)){
		if (sum(is.nan(allpower[,j])) == 0){
			if (is.vector(nobs)){
			plot(nobs, allpower[, j], ylab=paste('Power of', pnames[j]), xlab='Sample size')
			lines(nobs, allpower[, j])
			}else{
			N<-apply(nobs, 1, sum)
			plot(N, allpower[, j], ylab=paste('Power of', pnames[j]), xlab='Sample size')
			lines(N, allpower[, j])
			}
		}
	}
	invisible(allpower)
}