lab:projects:18non-normal_power_analysis:index
Some results
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))
R 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")
}
}