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))
## 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")
}
}