Nothing
R/bmem.R
In bmem: Mediation analysis with missing data using bootstrap
# Package: bmem
# Type: Package
# Title: Mediation analysis with missing data using bootstrap
# Version: 1.3
# Date: 2011-01-04
# Author: Zhiyong Zhang and Lijuan Wang
# Maintainer: Zhiyong Zhang <zhiyongzhang@nd.edu>
# Depends: R (>= 1.7), sem, Amelia, MASS
# Description: Four methods for mediation analysis with missing data: Listwise deletion, Pairwise deletion, Multiple imputation, and Two Stage Maximum Likelihood algorithm. For MI and TS-ML, auxiliary variables can be included. Bootstrap confidence intervals for mediation effects are obtained. The robust method is also implemented for TS-ML.
bmem.moments<-function(x, type=0){
#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, moment=FALSE){
## x: data
## moment=FALSE: covariance analysis only
listdata<-na.omit(x)
if (moment){
temp.cov<-bmem.moments(listdata, type=0)
}else{
temp.cov<-cov(listdata)
}
temp.cov
}
bmem.pair.cov<-function(x, moment=FALSE){
## x: data
## moment=FALSE: covariance analysis only
if (moment){
temp.cov<-bmem.moments(x, type=1)
}else{
temp.cov<-cov(x, use='pairwise.complete.obs')
}
temp.cov
}
bmem.sem<-function(x, ram, N, indirect, moment=FALSE, ...){
if (moment){
sem.res<-sem(ram, x, N, raw=TRUE, fixed.x='(i)', ...)
est<-sem.res$coeff
}else{
sem.res<-sem(ram, x, N, ...)
est<-sem.res$coeff
}
if (!missing(indirect)){
n<-length(indirect)
est.indirect<-rep(0,n)
for (i in 1:n){
temp<-indirect[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=indirect[i])
par.list<-as.list(est[par])
est.indirect[i]<-eval(ind.exp, par.list)
}
names(est.indirect)<-indirect
}else{
est.indirect<-NULL
}
model.fit<-NA
model.summary<-summary(sem.res)
model.fit<-c(model.summary$chisq,model.summary$GFI, model.summary$AGFI, model.summary$RMSEA[1], model.summary$NFI, model.summary$NNFI, model.summary$CFI, model.summary$BIC, model.summary$SRMR)
list(est=c(est, est.indirect), model.fit=model.fit)
}
bmem.sobel.ind<-function(sem.object, ind){
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, indirect, moment=FALSE, ...){
N<-nrow(x)
temp.cov<-bmem.list.cov(x, moment)
if (moment){
sem.object<-sem(ram, temp.cov, N, raw=TRUE, fixed.x='(i)', ...)
}else{
sem.object<-sem(ram, temp.cov, N, ...)
}
if (!missing(indirect)){
n<-length(indirect)
est.indirect<-NULL
for (i in 1:n){
temp<-indirect[i]
res<-bmem.sobel.ind(sem.object, temp)
est.indirect<-rbind(est.indirect, res)
}
sem.est<-summary(sem.object)$coeff[,1:3]
colnames(sem.est)<-colnames(est.indirect)
}else{
est.indirect<-NULL
sem.est<-summary(sem.object)$coeff[,1:3]
}
#print(summary(sem.object))
all.res<-rbind(sem.est, est.indirect)
p.value<-2*(1-pnorm(abs(all.res[,3])))
all.res<-cbind(all.res, p.value)
##cat('\n\nIndirect or mediation effects\n')
print(all.res)
invisible(list(estimates=all.res, sem=sem.object))
}
bmem.list<-function(x, ram, indirect, moment=FALSE, ...){
N<-dim(x)[1]
temp.cov<-bmem.list.cov(x, moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.pair<-function(x, ram, indirect, moment=FALSE, ...){
N<-dim(x)[1]
temp.cov<-bmem.pair.cov(x, moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.mi.cov<-function(x, m=10, moment=FALSE){
mi.data<-amelia(x,m,0)$imputations
mi.cov<-list()
for (j in 1:m){
temp.cov<-bmem.pair.cov(mi.data[[j]], moment)
mi.cov[[j]]<-temp.cov
}
mi.cov
}
bmem.v<-function(x, v, moment=FALSE){
if (moment){
v<-v+1
v<-c(1,v)
}
x[v,v]
}
bmem.mi<-function(x, ram, indirect, v, m=10, moment=FALSE, ...){
N<-dim(x)[1]
if (missing(v)) v<-1:(ncol(x))
temp.cov<-bmem.mi.cov(x,m,moment)
mi.est<-NULL
mi.fit<-NULL
for (j in 1:m){
v.cov<-bmem.v(temp.cov[[j]], v, moment)
temp<-bmem.sem(v.cov, ram, N, indirect, moment, ...)
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)
}
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, moment=FALSE, 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)
if (moment){
p1<-p+1
temp.cov<-matrix(NA, p+1,p+1)
temp.cov[1,1]<-1
temp.cov[1,2:p1]<-mu1
temp.cov[2:p1, 1]<-mu1
temp.cov[2:p1, 2:p1]<-sig1+m%x%t(m)
rownames(temp.cov)<-colnames(temp.cov)<-c('(i)', colnames(sig1))
}
temp.cov
}
bmem.em.rcov<-function(xmis, varphi=.1, moment=FALSE, 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)
if (moment){
p1<-p+1
temp.cov<-matrix(NA, p+1,p+1)
temp.cov[1,1]<-1
temp.cov[1,2:p1]<-mu1
temp.cov[2:p1, 1]<-mu1
temp.cov[2:p1, 2:p1]<-sig1+m%x%t(m)
rownames(temp.cov)<-colnames(temp.cov)<-c('(i)', colnames(x))
}
temp.cov
}
bmem.em<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', moment=FALSE, 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, moment, max_it, st)
}else{
temp.cov<-bmem.em.cov(temp.pattern, moment, max_it)
}
temp.cov<-bmem.v(temp.cov,v,moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.list.jack<-function(x, ram, indirect, moment=FALSE, ...){
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, indirect, moment, ...))
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)
}
bmem.pair.jack<-function(x, ram, indirect, moment=FALSE, ...){
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, indirect, moment, ...))
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)
}
bmem.mi.jack<-function(x, ram, indirect, v, m=10, moment=FALSE, ...){
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, indirect, v, m, moment, ...))
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)
}
bmem.em.jack<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', moment=FALSE, 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, indirect, v, robust, varphi, st, moment, 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)
}
bmem.list.boot<-function(x, ram, indirect, boot=1000, moment=FALSE, ...){
model0<-bmem.list(x, ram, indirect, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n,n, replace=TRUE),]
modelb<-try(bmem.list(x.boot, ram, indirect, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.pair.boot<-function(x, ram, indirect, boot=1000, moment=FALSE, ...){
model0<-bmem.pair(x, ram, indirect, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.pair(x.boot, ram, indirect, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.mi.boot<-function(x, ram, indirect, v, m=10, boot=1000, moment=FALSE, ...){
model0<-bmem.mi(x, ram, indirect, v, m, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.mi(x.boot, ram, indirect, v, m, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.em.boot<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', boot=1000, moment=FALSE, max_it=500, ...){
model0<-bmem.em(x, ram, indirect, v, robust, varphi, st, moment, max_it, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.em(x.boot, ram, indirect, v, robust, varphi, st, moment, max_it, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
## 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(x, ram, indirect, v, method='tsml', ci='bc', cl=.95, boot=1000, m=10, varphi=.1, st='i', robust=FALSE, max_it=500, moment=FALSE, ...){
## method: list, pair, mi, me
## ci: norm, perc, bc, bca
N<-nrow(x)
P<-ncol(x)
if (missing(v)) v<-1:P
## for listwise deletion method
if (method=='list'){
x<-x[,v]
boot.est<-bmem.list.boot(x, ram, indirect, boot, moment, ...)
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(x,ram,indirect,moment,...)
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'){
x<-x[,v]
boot.est<-bmem.pair.boot(x, ram, indirect, boot, moment, ...)
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(x,ram,indirect,moment,...)
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(x,ram,indirect,v,m,boot,moment,...)
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(x,ram,indirect,v,m,moment,...)
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(x,ram, indirect, v, robust, varphi, st, boot, moment, 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, indirect, v, robust, varphi, st, moment, 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')
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)
}
summary.bmem<-function(object, ci='bc', cl=.95, ...){
boot.est<-object$boot.est
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<-object$jack.est
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)
allci<-list(ci.est=ci.est, ci.fit=ci.fit)
class(allci)<-'summary.bmem'
invisible(allci)
}
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, indirect, v, ci='bc', cl=.95, boot=1000, max_it=500, ...){
## Estimate the saturated mean and covariance matrix
moment<-TRUE
xmiss<-bmem.pattern(x)
s.cov<-bmem.em.cov(xmiss, moment=TRUE, 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, ...)
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,indirect,v,moment,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,moment=FALSE, debug=FALSE){
if (moment){
fixed.x<-'(i)'
}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
}
## Monte Calor based power analysis for indirect effect using bootstrap
## by Zhiyong Zhang
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/Mediation effects:\n")
NAMES[def.idx] <- makeNames( object@ParTable$lhs[def.idx], "")
for(i in def.idx) {
print.estimate(name=NAMES[i], i)
}
cat("\n")
}
}
## bootstrap confidence intervals
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]<-1
}else{
par.value[i]<-0
}
}
}
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=100, nrep=1000, alpha=.95, skewness=NULL, kurtosis=NULL, ovnames=NULL, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
if (missing(model)) stop("A model is needed.")
model.indirect<-paste(model, "\n", indirect, "\n")
ngroups <- length(nobs)
## Initial analysis for some model information
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...)
if (ngroups > 1){
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', ...)
}else{
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
}
par.value<-popPar(temp.res)
idx <- 1:length( temp.res@ParTable$lhs )
#cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
ov<-lavNames(temp.res,'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]
}
newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
out<-temp.res
runonce<-function(i){
library(lavaan)
## Step 1: generate data
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
## Step 2: fit the model
if (ngroups > 1){
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...))
}else{
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...))
}
## Step 3: Check significance
if (class(temp.res)!="try-error"){
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])
}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)
have_mc <- have_snow <- FALSE
ncore <- ncore
if (parallel != "no" && ncore > 1L) {
if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
else if (parallel == "snow") have_snow <- TRUE
if (!have_mc && !have_snow) ncore <- 1L
}
RR <- nrep
res <- if (ncore > 1L && (have_mc || have_snow)) {
if (have_mc) {
parallel::mclapply(seq_len(RR), runonce, mc.cores = ncore)
} else if (have_snow) {
list(...) # evaluate any promises
if (is.null(cl)) {
cl <- parallel::makePSOCKcluster(rep("localhost", ncore))
if(RNGkind()[1L] == "L'Ecuyer-CMRG")
parallel::clusterSetRNGStream(cl)
res <- parallel::parLapply(cl, seq_len(RR), runonce)
parallel::stopCluster(cl)
res
} else parallel::parLapply(cl, 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)
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, 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, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
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
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...)
if (ngroups > 1){
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...)
}else{
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...)
}
par.value<-popPar(temp.res)
idx <- 1:length( temp.res@ParTable$lhs )
#cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
ov<-lavNames(temp.res,'ov')
ptype<-temp.res@ParTable$op
## 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]
}
newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
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)
}
## repeated the following for nrep times
runonce<-function(i){
library(lavaan)
## Step 1: generate data
newdata<-lavaan:::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, ...)
if (ngroups > 1){
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...))
}else{
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...))
}
## Step 3: Conduct bootstrap analysis
if (class(temp.res)!="try-error"){
orig.res<-coef.new(temp.res)
boot.res<-lavaan:::bootstrapLavaan(temp.res, FUN=coef.new, R=nboot, parallel="no", warn=FALSE, ...)
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
# this is from the boot function in package boot
old_options <- options(); options(warn = -1)
have_mc <- have_snow <- FALSE
ncore <- ncore
if (parallel != "no" && ncore > 1L) {
if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
else if (parallel == "snow") have_snow <- TRUE
if (!have_mc && !have_snow) ncore <- 1L
}
RR <- nrep
res <- if (ncore > 1L && (have_mc || have_snow)) {
if (have_mc) {
parallel::mclapply(seq_len(RR), runonce, mc.cores = ncore)
} else if (have_snow) {
list(...) # evaluate any promises
if (is.null(cl)) {
cl <- parallel::makePSOCKcluster(rep("localhost", ncore))
if(RNGkind()[1L] == "L'Ecuyer-CMRG")
parallel::clusterSetRNGStream(cl)
res <- parallel::parLapply(cl, seq_len(RR), runonce)
parallel::stopCluster(cl)
res
} else parallel::parLapply(cl, 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)
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, 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, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
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, cl=cl, ...)
}else{
indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore, cl=cl, ...)
}
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, cl=cl, ...)
}else{
indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore, cl=cl, ...)
}
allpower<-rbind(allpower, indpower$power)
}
}
windows(record=TRUE)
op <- par(ask=TRUE)
on.exit(par(op))
pnames <- colnames(allpower)
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)
}
Try the bmem package in your browser
Any scripts or data that you put into this service are public.
bmem documentation built on May 2, 2019, 5:27 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Package: bmem
# Type: Package
# Title: Mediation analysis with missing data using bootstrap
# Version: 1.3
# Date: 2011-01-04
# Author: Zhiyong Zhang and Lijuan Wang
# Maintainer: Zhiyong Zhang <zhiyongzhang@nd.edu>
# Depends: R (>= 1.7), sem, Amelia, MASS
# Description: Four methods for mediation analysis with missing data: Listwise deletion, Pairwise deletion, Multiple imputation, and Two Stage Maximum Likelihood algorithm. For MI and TS-ML, auxiliary variables can be included. Bootstrap confidence intervals for mediation effects are obtained. The robust method is also implemented for TS-ML.
bmem.moments<-function(x, type=0){
#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, moment=FALSE){
## x: data
## moment=FALSE: covariance analysis only
listdata<-na.omit(x)
if (moment){
temp.cov<-bmem.moments(listdata, type=0)
}else{
temp.cov<-cov(listdata)
}
temp.cov
}
bmem.pair.cov<-function(x, moment=FALSE){
## x: data
## moment=FALSE: covariance analysis only
if (moment){
temp.cov<-bmem.moments(x, type=1)
}else{
temp.cov<-cov(x, use='pairwise.complete.obs')
}
temp.cov
}
bmem.sem<-function(x, ram, N, indirect, moment=FALSE, ...){
if (moment){
sem.res<-sem(ram, x, N, raw=TRUE, fixed.x='(i)', ...)
est<-sem.res$coeff
}else{
sem.res<-sem(ram, x, N, ...)
est<-sem.res$coeff
}
if (!missing(indirect)){
n<-length(indirect)
est.indirect<-rep(0,n)
for (i in 1:n){
temp<-indirect[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=indirect[i])
par.list<-as.list(est[par])
est.indirect[i]<-eval(ind.exp, par.list)
}
names(est.indirect)<-indirect
}else{
est.indirect<-NULL
}
model.fit<-NA
model.summary<-summary(sem.res)
model.fit<-c(model.summary$chisq,model.summary$GFI, model.summary$AGFI, model.summary$RMSEA[1], model.summary$NFI, model.summary$NNFI, model.summary$CFI, model.summary$BIC, model.summary$SRMR)
list(est=c(est, est.indirect), model.fit=model.fit)
}
bmem.sobel.ind<-function(sem.object, ind){
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, indirect, moment=FALSE, ...){
N<-nrow(x)
temp.cov<-bmem.list.cov(x, moment)
if (moment){
sem.object<-sem(ram, temp.cov, N, raw=TRUE, fixed.x='(i)', ...)
}else{
sem.object<-sem(ram, temp.cov, N, ...)
}
if (!missing(indirect)){
n<-length(indirect)
est.indirect<-NULL
for (i in 1:n){
temp<-indirect[i]
res<-bmem.sobel.ind(sem.object, temp)
est.indirect<-rbind(est.indirect, res)
}
sem.est<-summary(sem.object)$coeff[,1:3]
colnames(sem.est)<-colnames(est.indirect)
}else{
est.indirect<-NULL
sem.est<-summary(sem.object)$coeff[,1:3]
}
#print(summary(sem.object))
all.res<-rbind(sem.est, est.indirect)
p.value<-2*(1-pnorm(abs(all.res[,3])))
all.res<-cbind(all.res, p.value)
##cat('\n\nIndirect or mediation effects\n')
print(all.res)
invisible(list(estimates=all.res, sem=sem.object))
}
bmem.list<-function(x, ram, indirect, moment=FALSE, ...){
N<-dim(x)[1]
temp.cov<-bmem.list.cov(x, moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.pair<-function(x, ram, indirect, moment=FALSE, ...){
N<-dim(x)[1]
temp.cov<-bmem.pair.cov(x, moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.mi.cov<-function(x, m=10, moment=FALSE){
mi.data<-amelia(x,m,0)$imputations
mi.cov<-list()
for (j in 1:m){
temp.cov<-bmem.pair.cov(mi.data[[j]], moment)
mi.cov[[j]]<-temp.cov
}
mi.cov
}
bmem.v<-function(x, v, moment=FALSE){
if (moment){
v<-v+1
v<-c(1,v)
}
x[v,v]
}
bmem.mi<-function(x, ram, indirect, v, m=10, moment=FALSE, ...){
N<-dim(x)[1]
if (missing(v)) v<-1:(ncol(x))
temp.cov<-bmem.mi.cov(x,m,moment)
mi.est<-NULL
mi.fit<-NULL
for (j in 1:m){
v.cov<-bmem.v(temp.cov[[j]], v, moment)
temp<-bmem.sem(v.cov, ram, N, indirect, moment, ...)
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)
}
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, moment=FALSE, 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)
if (moment){
p1<-p+1
temp.cov<-matrix(NA, p+1,p+1)
temp.cov[1,1]<-1
temp.cov[1,2:p1]<-mu1
temp.cov[2:p1, 1]<-mu1
temp.cov[2:p1, 2:p1]<-sig1+m%x%t(m)
rownames(temp.cov)<-colnames(temp.cov)<-c('(i)', colnames(sig1))
}
temp.cov
}
bmem.em.rcov<-function(xmis, varphi=.1, moment=FALSE, 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)
if (moment){
p1<-p+1
temp.cov<-matrix(NA, p+1,p+1)
temp.cov[1,1]<-1
temp.cov[1,2:p1]<-mu1
temp.cov[2:p1, 1]<-mu1
temp.cov[2:p1, 2:p1]<-sig1+m%x%t(m)
rownames(temp.cov)<-colnames(temp.cov)<-c('(i)', colnames(x))
}
temp.cov
}
bmem.em<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', moment=FALSE, 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, moment, max_it, st)
}else{
temp.cov<-bmem.em.cov(temp.pattern, moment, max_it)
}
temp.cov<-bmem.v(temp.cov,v,moment)
bmem.sem(temp.cov, ram, N, indirect, moment, ...)
}
bmem.list.jack<-function(x, ram, indirect, moment=FALSE, ...){
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, indirect, moment, ...))
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)
}
bmem.pair.jack<-function(x, ram, indirect, moment=FALSE, ...){
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, indirect, moment, ...))
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)
}
bmem.mi.jack<-function(x, ram, indirect, v, m=10, moment=FALSE, ...){
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, indirect, v, m, moment, ...))
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)
}
bmem.em.jack<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', moment=FALSE, 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, indirect, v, robust, varphi, st, moment, 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)
}
bmem.list.boot<-function(x, ram, indirect, boot=1000, moment=FALSE, ...){
model0<-bmem.list(x, ram, indirect, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n,n, replace=TRUE),]
modelb<-try(bmem.list(x.boot, ram, indirect, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.pair.boot<-function(x, ram, indirect, boot=1000, moment=FALSE, ...){
model0<-bmem.pair(x, ram, indirect, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.pair(x.boot, ram, indirect, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.mi.boot<-function(x, ram, indirect, v, m=10, boot=1000, moment=FALSE, ...){
model0<-bmem.mi(x, ram, indirect, v, m, moment, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.mi(x.boot, ram, indirect, v, m, moment, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
bmem.em.boot<-function(x, ram, indirect, v, robust=FALSE, varphi=.1, st='i', boot=1000, moment=FALSE, max_it=500, ...){
model0<-bmem.em(x, ram, indirect, v, robust, varphi, st, moment, max_it, ...)
par0<-model0$est
fit0<-model0$model.fit
n<-dim(x)[1]
boot.est<-NULL
boot.fit<-NULL
for (i in 1:boot){
x.boot<-x[sample(n, n, replace=TRUE),]
modelb<-try(bmem.em(x.boot, ram, indirect, v, robust, varphi, st, moment, max_it, ...))
if (class(modelb)!="try-error"){
boot.est<-rbind(boot.est, modelb$est)
boot.fit<-rbind(boot.fit, modelb$model.fit)
}
}
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)
}
## 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(x, ram, indirect, v, method='tsml', ci='bc', cl=.95, boot=1000, m=10, varphi=.1, st='i', robust=FALSE, max_it=500, moment=FALSE, ...){
## method: list, pair, mi, me
## ci: norm, perc, bc, bca
N<-nrow(x)
P<-ncol(x)
if (missing(v)) v<-1:P
## for listwise deletion method
if (method=='list'){
x<-x[,v]
boot.est<-bmem.list.boot(x, ram, indirect, boot, moment, ...)
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(x,ram,indirect,moment,...)
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'){
x<-x[,v]
boot.est<-bmem.pair.boot(x, ram, indirect, boot, moment, ...)
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(x,ram,indirect,moment,...)
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(x,ram,indirect,v,m,boot,moment,...)
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(x,ram,indirect,v,m,moment,...)
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(x,ram, indirect, v, robust, varphi, st, boot, moment, 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, indirect, v, robust, varphi, st, moment, 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')
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)
}
summary.bmem<-function(object, ci='bc', cl=.95, ...){
boot.est<-object$boot.est
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<-object$jack.est
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)
allci<-list(ci.est=ci.est, ci.fit=ci.fit)
class(allci)<-'summary.bmem'
invisible(allci)
}
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, indirect, v, ci='bc', cl=.95, boot=1000, max_it=500, ...){
## Estimate the saturated mean and covariance matrix
moment<-TRUE
xmiss<-bmem.pattern(x)
s.cov<-bmem.em.cov(xmiss, moment=TRUE, 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, ...)
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,indirect,v,moment,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,moment=FALSE, debug=FALSE){
if (moment){
fixed.x<-'(i)'
}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
}
## Monte Calor based power analysis for indirect effect using bootstrap
## by Zhiyong Zhang
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/Mediation effects:\n")
NAMES[def.idx] <- makeNames( object@ParTable$lhs[def.idx], "")
for(i in def.idx) {
print.estimate(name=NAMES[i], i)
}
cat("\n")
}
}
## bootstrap confidence intervals
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]<-1
}else{
par.value[i]<-0
}
}
}
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=100, nrep=1000, alpha=.95, skewness=NULL, kurtosis=NULL, ovnames=NULL, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
if (missing(model)) stop("A model is needed.")
model.indirect<-paste(model, "\n", indirect, "\n")
ngroups <- length(nobs)
## Initial analysis for some model information
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...)
if (ngroups > 1){
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', ...)
}else{
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, ...)
}
par.value<-popPar(temp.res)
idx <- 1:length( temp.res@ParTable$lhs )
#cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
ov<-lavNames(temp.res,'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]
}
newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
out<-temp.res
runonce<-function(i){
library(lavaan)
## Step 1: generate data
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
## Step 2: fit the model
if (ngroups > 1){
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...))
}else{
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...))
}
## Step 3: Check significance
if (class(temp.res)!="try-error"){
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])
}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)
have_mc <- have_snow <- FALSE
ncore <- ncore
if (parallel != "no" && ncore > 1L) {
if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
else if (parallel == "snow") have_snow <- TRUE
if (!have_mc && !have_snow) ncore <- 1L
}
RR <- nrep
res <- if (ncore > 1L && (have_mc || have_snow)) {
if (have_mc) {
parallel::mclapply(seq_len(RR), runonce, mc.cores = ncore)
} else if (have_snow) {
list(...) # evaluate any promises
if (is.null(cl)) {
cl <- parallel::makePSOCKcluster(rep("localhost", ncore))
if(RNGkind()[1L] == "L'Ecuyer-CMRG")
parallel::clusterSetRNGStream(cl)
res <- parallel::parLapply(cl, seq_len(RR), runonce)
parallel::stopCluster(cl)
res
} else parallel::parLapply(cl, 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)
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, 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, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
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
newdata<-lavaan:::simulateData(model,sample.nobs=nobs,skewness=0,kurtosis=0, ...)
if (ngroups > 1){
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...)
}else{
temp.res<-lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...)
}
par.value<-popPar(temp.res)
idx <- 1:length( temp.res@ParTable$lhs )
#cnames<-paste(temp.res@ParTable$lhs[idx], temp.res@ParTable$op[idx], temp.res@ParTable$rhs[idx])
ov<-lavNames(temp.res,'ov')
ptype<-temp.res@ParTable$op
## 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]
}
newdata<-simulateData(model,sample.nobs=nobs,skewness=skewness,kurtosis=kurtosis, ...)
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)
}
## repeated the following for nrep times
runonce<-function(i){
library(lavaan)
## Step 1: generate data
newdata<-lavaan:::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, ...)
if (ngroups > 1){
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, group='group', warn=FALSE, ...))
}else{
temp.res<-try(lavaan:::sem(model.indirect, data=newdata, se=se, estimator=estimator, warn=FALSE, ...))
}
## Step 3: Conduct bootstrap analysis
if (class(temp.res)!="try-error"){
orig.res<-coef.new(temp.res)
boot.res<-lavaan:::bootstrapLavaan(temp.res, FUN=coef.new, R=nboot, parallel="no", warn=FALSE, ...)
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
# this is from the boot function in package boot
old_options <- options(); options(warn = -1)
have_mc <- have_snow <- FALSE
ncore <- ncore
if (parallel != "no" && ncore > 1L) {
if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
else if (parallel == "snow") have_snow <- TRUE
if (!have_mc && !have_snow) ncore <- 1L
}
RR <- nrep
res <- if (ncore > 1L && (have_mc || have_snow)) {
if (have_mc) {
parallel::mclapply(seq_len(RR), runonce, mc.cores = ncore)
} else if (have_snow) {
list(...) # evaluate any promises
if (is.null(cl)) {
cl <- parallel::makePSOCKcluster(rep("localhost", ncore))
if(RNGkind()[1L] == "L'Ecuyer-CMRG")
parallel::clusterSetRNGStream(cl)
res <- parallel::parLapply(cl, seq_len(RR), runonce)
parallel::stopCluster(cl)
res
} else parallel::parLapply(cl, 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)
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, 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, se="default", estimator="default", parallel="no", ncore=Sys.getenv('NUMBER_OF_PROCESSORS'), cl=NULL, ...){
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, cl=cl, ...)
}else{
indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore, cl=cl, ...)
}
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, cl=cl, ...)
}else{
indpower <- power.boot(model=model, indirect=indirect, nobs=N, nrep=nrep, nboot=nboot, alpha=alpha, skewness=skewness, kurtosis=kurtosis, ovnames=ovnames, se=se, estimator=estimator, parallel=parallel, ncore=ncore, cl=cl, ...)
}
allpower<-rbind(allpower, indpower$power)
}
}
windows(record=TRUE)
op <- par(ask=TRUE)
on.exit(par(op))
pnames <- colnames(allpower)
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)
}
Try the bmem package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.