Nothing
## The package call the rsem package for robust analysis
## R function for Cronbach alpha
alpha<-function(y, varphi=0.1, se=FALSE, test=TRUE, complete=FALSE, auxiliary=NULL, drop, silent=TRUE){
## get the dim of y
p<-ncol(y)
## combine data
if (!is.null(auxiliary)){
if (!(is.matrix(auxiliary) | is.data.frame(auxiliary))) stop('auxiliary should be an matrix or a data frame')
y.aux<-cbind(y, auxiliary)
p.aux<-ncol(auxiliary)
}else{
y.aux<-y
}
## remove some cases
if (!missing(drop)){
drop <- -drop
obsind <- 1:nrow(y.aux)
obsind <- obsind[drop]
y.aux<-y.aux[obsind, ]
}
if (complete) y.aux<-na.omit(y.aux)
n<-nrow(y.aux)
rownames(y.aux)<-1:n
pat<-rsem.pattern(y.aux)
cov1<-rsem.emmusig(pat,varphi=varphi)
weight<-rep(1, n)
if (varphi>0) weight<-cov1$weight
## covariance for y only
sigma<-cov1$sigma[1:p, 1:p]
mu<-cov1$mu[1:p]
alpha<-p/(p-1)*(1-sum(diag(sigma))/sum(sigma))
## robust standard error
se.alpha<-NA
if (se | test){
s<-rsem.Ascov(pat, cov1, varphi=varphi)
## get out the Gamma matrix
vname<-rownames(sigma)
gname<-NULL
for (i in 1:p){
for (j in i:p){
gname<-c(gname, paste(vname[i], ".", vname[j], sep=""))
}
}
gmname<-c(vname, gname)
}
if (se){
firstd<-NULL
for (i in 1:p){
for (j in i:p){
if (i==j){
temp<--p/(p-1)*(sum(sigma)-sum(diag(sigma)))/((sum(sigma))^2)
}else{
temp<- 2*p/(p-1)*sum(diag(sigma))/((sum(sigma))^2)
}
firstd<-c(firstd, temp)
}
}
## take out the part of gamma for the covariance of y
if (is.null(auxiliary)){
q<-p+p*(p+1)/2
gamma<-s$Gamma[(p+1):q, (p+1):q]
}else{
gamma<-s$Gamma[gname, gname]
}
se.alpha<-sqrt(firstd%*%gamma%*%t(t(firstd)))/sqrt(n)
}
## conduct the test of tau-equivalent
fit<-NULL
if (test){
if (is.null(auxiliary)){
gamma<-s$Gamma
}else{
gamma<-s$Gamma[gmname, gmname]
}
vname<-rownames(sigma)
## model for tau-equivalent
model<-paste('f =~ 1*', vname[1])
for (i in 2:length(vname)){
model<-paste(model, '+ 1*', vname[i])
}
## estimate the model using rsem
cfa.res<-cfa(model, sample.cov=sigma, sample.mean=mu, meanstructure = TRUE, sample.nobs=n)
robust.fit <- rsem.fit(cfa.res, gamma, cov1)
fit<-robust.fit[[4]][[1]]
if (silent){
cat("Test of tau-equavilence\n")
cat("The robust F statistic is ", round(fit[1],3), "\n")
cat("with a p-value ", round(fit[4],4), "\n")
if (fit[4]<.05) cat("**The F test rejected tau-equavilence**\n\n")
}
}
if (varphi>0) prop.down<-(n-sum(weight$w1==1))/n*100
if (silent) cat('The alpha is ', alpha, sep='')
if (se & silent) cat(' with the standard error ', se.alpha, sep='')
if (silent) cat('.\n')
if (varphi>0 & silent) cat('About ', round(prop.down,2), '% of cases were downweighted.\n', sep='')
alphaobject<-list(alpha=alpha, se=se.alpha, weight=weight, y=y, auxiliary=auxiliary, varphi=varphi, musig=cov1, fit=fit)
class(alphaobject)<-'alpha'
invisible(alphaobject)
}
plot.alpha<-function(x, type="weight", profile=5, interval=0.01, center=TRUE, scale=FALSE, w1=FALSE, numbered=FALSE, pos='topright', ...){
## type: weight, profile, diagnosis
res<-x
y<-res$y
outcase.temp<-sum(res$weight$w1<1)
## find the smallest weight
if (profile==0){
outcase<-5
if (outcase.temp<5) outcase<-outcase.temp
}
if (outcase.temp<5) profile<-outcase.temp
## find the smallest cases
temp<-sort(res$weight$w1)
if (profile==0){
idx<-which(res$weight$w1 <= temp[outcase])
}else{
idx<-which(res$weight$w1 <= temp[profile])
}
if (substr(type,1,1)=="w"){
if (w1){
par(mfrow=c(2,1))
plot(res$weight$w1, ylim=c(0,1), xlab='Case number', ylab='Weight w1', main='Weights for mean estimates',...)
text(idx, res$weight$w1[idx], idx, pos=1, col='red', ...)
plot(res$weight$w2, ylim=c(0,max(res$weight$w2)), xlab='Case number', ylab='Weight w2', main='Weights for covariance estimates',...)
text(idx, res$weight$w2[idx], idx, pos=1, col='red', ...)
par(mfrow=c(1,1))
}else{
plot(res$weight$w2, ylim=c(0,max(res$weight$w2)), xlab='Case number', ylab='Weight w2',...)
text(idx, res$weight$w2[idx], idx, pos=1, col='red', ...)
}
}
if (substr(type,1,1)=="p"){
## profile plot
if (center){
if (scale){
for (i in 1:nrow(y)){
y[i, ]<-(y[i, ]-res$musig$mu)/sqrt(diag(res$musig$sigma))
}
}else{
for (i in 1:nrow(y)){
y[i, ]<-y[i, ]-res$musig$mu
}
}
}
l.min<-min(y, na.rm=TRUE)
l.max<-max(y, na.rm=TRUE)
p<-ncol(y)
if (center){
plot(1:p, rep(0,p), type='l', ylim=c(l.min, l.max), lwd=3, ylab='Score', xlab='', axes=FALSE, col='black', ...)
}else{
plot(1:p, res$musig$mu, type='l', ylim=c(l.min, l.max), lwd=3, ylab='Score', xlab='', axes=FALSE, col='black', ...)
}
box()
axis(1, at=1:p, labels=names(y), las=2, ...)
axis(2, ...)
if (!is.null(pos)){
ltyno<-1
idx.type<-NULL
for (i in idx){
lines(1:p, y[i, ], lty=ltyno, ...)
if (numbered & center) text(1, y[i,1], i)
ltyno<-ltyno+1
## type of outliers
temp.type<-'(O)'
if (max(y[i, ], na.rm=TRUE)<0) temp.type<-'(L-)'
if (min(y[i, ], na.rm=TRUE)>0) temp.type<-'(L+)'
idx.type<- c(idx.type, temp.type)
}
legend(pos, legend=paste(idx, idx.type), lty=1:ltyno)
}
}
if (substr(type,1,1)=="d"){
varphi<-res$varphi
phi<-seq(0, varphi, by=interval)
alpha.diag<-NULL
for (i in 1:length(phi)){
alpha.diag<-c(alpha.diag, alpha(y, phi[i], auxiliary=res$auxiliary, test=FALSE, se=FALSE)$alpha)
}
plot(phi, alpha.diag, type='o', xlab='varphi', ylab='alpha', ylim=c(0,1), ...)
points(phi, alpha.diag, bg='black', pch=21, cex=1, lwd=1)
}
}
summary.alpha<-function(object, type='raw', prob=.95, ...){
if (prob > .5) prob <- 1 - prob
res<-object
cat("The estimated alpha is \n")
t0.txt <- sprintf(" %-20s", "alpha")
t1.txt <- sprintf(" %10.3f", res$alpha)
cat(t0.txt, t1.txt, "\n", sep="")
if (!is.na(res$se)){
t0.txt <- sprintf(" %-20s", "se")
t1.txt <- sprintf(" %10.3f", res$se)
cat(t0.txt, t1.txt, "\n", sep="")
## output p-value, too
z<-abs(res$alpha/res$se)
pvalue<-(1-pnorm(z))
t0.txt <- sprintf(" %-20s", "p-value (alpha>0)")
t1.txt <- sprintf(" %10.3f", pvalue)
cat(t0.txt, t1.txt, "\n", sep="")
## output confidence interval
if (type=='logit'){
lambda<-log(res$alpha/(1-res$alpha))
se.lambda<-res$se/(res$alpha*(1-res$alpha))
ci<-c( 1/(1+exp(-lambda+abs(qnorm(prob/2))*se.lambda)), 1/(1+exp(-lambda-abs(qnorm(prob/2))*se.lambda)))
}else{
ci<-c(res$alpha-abs(qnorm(prob/2))*res$se, res$alpha+abs(qnorm(prob/2))*res$se)
if (ci[2]>1) ci[2]<-1.000
if (ci[1]<0) ci[1]<-0.000
}
t0.txt <- sprintf(" %-20s", "Confidence interval")
t1.txt <- sprintf(" %10.3f", ci)
cat(t0.txt, t1.txt, "\n\n", sep="")
}
## output the fit statistic if existing
if (!is.null(res$fit)){
cat("\nTest of tau-equivalence \n")
cat(" The robust F statistic is ", round(res$fit[1],3), "\n")
cat(" with a p-value ", round(res$fit[4],4), "\n")
if (res$fit[4]<.05){
cat(" Note. The robust test rejected the tau-equivalence assumption.\n")
}else{
cat(" Note. The robust test failed to reject the tau-equivalence assumption.\n")
}
}
}
omega<-function(y, varphi=0.1, se=FALSE, test=TRUE, complete=FALSE, auxiliary=NULL, drop, silent=TRUE){
## get the dim of y
p<-ncol(y)
## combine data
if (!is.null(auxiliary)){
if (!(is.matrix(auxiliary) | is.data.frame(auxiliary))) stop('auxiliary should be an matrix or a data frame')
y.aux<-cbind(y, auxiliary)
p.aux<-ncol(auxiliary)
}else{
y.aux<-y
}
## remove some cases
if (!missing(drop)){
drop <- -drop
obsind <- 1:nrow(y.aux)
obsind <- obsind[drop]
y.aux<-y.aux[obsind, ]
}
if (complete) y.aux<-na.omit(y.aux)
n<-nrow(y.aux)
rownames(y.aux)<-1:n
pat<-rsem.pattern(y.aux)
cov1<-rsem.emmusig(pat,varphi=varphi)
weight<-rep(1, n)
if (varphi>0) weight<-cov1$weight
## covariance for y only
sigma<-cov1$sigma[1:p, 1:p]
mu<-cov1$mu[1:p]
## fit a factor model with the covariance matrix and get parameter estimates and their covariance matrix
vname<-rownames(sigma)
## the factor model
model<-paste('f =~ ', vname[1])
for (i in 2:length(vname)){
model<-paste(model, '+', vname[i])
}
## estimate the model using rsem
cfa.res<-cfa(model, sample.cov=sigma, sample.mean=cov1$mu, meanstructure = TRUE, sample.nobs=n, std.lv=TRUE)
## calculate omega
partable <- parTable(cfa.res)
parest <- parameterEstimates(cfa.res)
cfa.load<-parest$est[partable$op=="=~" & partable$free>0]
cfa.psi<-parest$est[partable$op=="~~" & partable$free>0]
p1<-(sum(cfa.load))^2
p2<- p1 + sum(cfa.psi)
omega<- p1/p2
## robust standard error
se.omega<-NA
if (se | test){
s<-rsem.Ascov(pat, cov1, varphi=varphi)
## get out the Gamma matrix
vname<-rownames(sigma)
gname<-NULL
for (i in 1:p){
for (j in i:p){
gname<-c(gname, paste(vname[i], ".", vname[j], sep=""))
}
}
gmname<-c(vname, gname)
Gamma<-s$Gamma[gmname, gmname]
}
if (se){
cfa.se<-rsem.se(cfa.res, Gamma)
firstd<-rep(sum(cfa.load)*sum(cfa.psi)*2/(p2^2) , p)
firstd<-c(firstd, rep( -p1/(p2^2),p))
q<-p*2
gamma<-cfa.se$vcov[[1]][1:q, 1:q]
se.omega<-sqrt(firstd%*%gamma%*%t(t(firstd)))/sqrt(n)
}
fit<-NULL
if (test){
robust.fit <- rsem.fit(cfa.res, Gamma, cov1)
fit<-robust.fit[[4]][[1]]
if (silent){
cat("Test of homogeneity\n")
cat("The robust F statistic is ", round(fit[1],3), "\n")
cat("with a p-value ", round(fit[4],4), "\n\n")
if (fit[4]<.05) cat("**The F test rejected homogeneity**\n\n")
}
}
if (varphi>0) prop.down<-(n-sum(weight$w1==1))/n*100
if (silent)cat('The omega is ', omega, sep='')
if (se & silent) cat(' with the standard error ', se.omega, sep='')
if (silent)cat('.\n')
if (varphi>0 & silent) cat('About ', round(prop.down,2), '% of cases were downweighted.\n', sep='')
omegaobject<-list(omega=omega, se=se.omega, weight=weight, y=y, auxiliary=auxiliary, varphi=varphi, musig=cov1, fit=fit)
class(omegaobject)<-'omega'
invisible(omegaobject)
}
summary.omega<-function(object, type='raw', prob=.95,...){
if (prob > .5) prob <- 1 - prob
res<-object
cat("The estimated omega is \n")
t0.txt <- sprintf(" %-20s", "omega")
t1.txt <- sprintf(" %10.3f", res$omega)
cat(t0.txt, t1.txt, "\n", sep="")
if (!is.na(res$se)){
t0.txt <- sprintf(" %-20s", "se")
t1.txt <- sprintf(" %10.3f", res$se)
cat(t0.txt, t1.txt, "\n", sep="")
## output p-value, too
z<-abs(res$omega/res$se)
pvalue<-(1-pnorm(z))
t0.txt <- sprintf(" %-20s", "p-value (omega>0)")
t1.txt <- sprintf(" %10.3f", pvalue)
cat(t0.txt, t1.txt, "\n", sep="")
## output confidence interval
if (type=='logit'){
lambda<-log(res$omega/(1-res$omega))
se.lambda<-res$se/(res$omega*(1-res$omega))
ci<-c( 1/(1+exp(-lambda+abs(qnorm(prob/2))*se.lambda)), 1/(1+exp(-lambda-abs(qnorm(prob/2))*se.lambda)))
}else{
ci<-c(res$omega-abs(qnorm(prob/2))*res$se, res$omega+abs(qnorm(prob/2))*res$se)
if (ci[2]>1) ci[2]<-1.000
if (ci[1]<0) ci[1]<-0.000
}
t0.txt <- sprintf(" %-20s", "Confidence interval")
t1.txt <- sprintf(" %10.3f", ci)
cat(t0.txt, t1.txt, "\n\n", sep="")
}
if (!is.null(res$fit)){
cat("\nTest of homogeneity\n")
cat(" The robust F statistic is ", round(res$fit[1],3), "\n")
cat(" with a p-value ", round(res$fit[4],4), "\n")
if (res$fit[4]<.05){
cat(" Note. The robust test rejected the assumption of homogeneity.\n")
}else{
cat(" Note. The robust test failed to reject the assumption of homogeneity.\n")
}
}
}
plot.omega<-function(x, type="weight", profile=5, interval=0.01, center=TRUE, scale=FALSE, w1=FALSE, numbered=FALSE, pos='topright', ...){
## type: weight, profile, diagnosis
res<-x
y<-res$y
outcase.temp<-sum(res$weight$w1<1)
## find the smallest weight
if (profile==5){
outcase<-5
if (outcase.temp<5) outcase<-outcase.temp
}
if (outcase.temp<5) profile<-outcase.temp
## find the smallest cases
temp<-sort(res$weight$w1)
if (profile==5){
idx<-which(res$weight$w1 <= temp[outcase])
}else{
idx<-which(res$weight$w1 <= temp[profile])
}
if (substr(type,1,1)=="w"){
if (w1){
par(mfrow=c(2,1))
plot(res$weight$w1, ylim=c(0,1), xlab='Case number', ylab='Weight w1', main='Weights for mean estimates',...)
text(idx, res$weight$w1[idx], idx, pos=1, col='red', ...)
plot(res$weight$w2, ylim=c(0,max(res$weight$w2)), xlab='Case number', ylab='Weight w2', main='Weights for covariance estimates',...)
text(idx, res$weight$w2[idx], idx, pos=1, col='red', ...)
par(mfrow=c(1,1))
}else{
plot(res$weight$w2, ylim=c(0,max(res$weight$w2)), xlab='Case number', ylab='Weight w2',...)
text(idx, res$weight$w2[idx], idx, pos=1, col='red', ...)
}
}
if (substr(type,1,1)=="p"){
## profile plot
if (center){
if (scale){
for (i in 1:nrow(y)){
y[i, ]<-(y[i, ]-res$musig$mu)/sqrt(diag(res$musig$sigma))
}
}else{
for (i in 1:nrow(y)){
y[i, ]<-y[i, ]-res$musig$mu
}
}
}
l.min<-min(y, na.rm=TRUE)
l.max<-max(y, na.rm=TRUE)
p<-ncol(y)
if (center){
plot(1:p, rep(0,p), type='l', ylim=c(l.min, l.max), lwd=3, ylab='Score', xlab='', axes=FALSE, col='black', ...)
}else{
plot(1:p, res$musig$mu, type='l', ylim=c(l.min, l.max), lwd=3, ylab='Score', xlab='', axes=FALSE, col='black', ...)
}
box()
axis(1, at=1:p, labels=names(y), las=2, ...)
axis(2, ...)
if (!is.null(pos)){
ltyno<-1
idx.type<-NULL
for (i in idx){
lines(1:p, y[i, ], lty=ltyno, ...)
if (numbered & center) text(1, y[i,1], i)
ltyno<-ltyno+1
## type of outliers
temp.type<-'(O)'
if (max(y[i, ], na.rm=TRUE)<0) temp.type<-'(L-)'
if (min(y[i, ], na.rm=TRUE)>0) temp.type<-'(L+)'
idx.type<- c(idx.type, temp.type)
}
legend(pos, legend=paste(idx, idx.type), lty=1:ltyno)
}
}
if (substr(type,1,1)=="d"){
varphi<-res$varphi
phi<-seq(0, varphi, by=interval)
omega.diag<-NULL
for (i in 1:length(phi)){
omega.diag<-c(omega.diag, omega(y, phi[i], auxiliary=res$auxiliary, test=FALSE, se=FALSE)$omega)
}
plot(phi, omega.diag, xlab='varphi', type='o', lwd=1, ylab='omega', ylim=c(0,1), ...)
points(phi, omega.diag, bg='black', pch=21, cex=1, lwd=1)
}
}
tau.test<-function(y, varphi=0.1, complete=FALSE, drop){
if (!missing(drop)){
drop <- -drop
obsind <- 1:nrow(y)
obsind <- obsind[drop]
y<-y[obsind, ]
}
if (complete) y<-na.omit(y)
p<-ncol(y)
n<-nrow(y)
rownames(y)<-1:n
pat<-rsem.pattern(y)
cov1<-rsem.emmusig(pat,varphi=varphi)
ascov<-rsem.Ascov(pat, cov1, varphi=varphi)
sigma<-cov1$sigma
## fit a factor model with the covariance matrix and get parameter estimates and their covariance matrix
vname<-rownames(sigma)
## estimate two models
## model for tau-equivalent
model<-paste('f =~ 1*', vname[1])
for (i in 2:length(vname)){
model<-paste(model, '+ 1*', vname[i])
}
## estimate the model using rsem
cfa.res1<-cfa(model, sample.cov=sigma, sample.mean=cov1$mu, meanstructure = TRUE, sample.nobs=n)
robust.fit1 <- rsem.fit(cfa.res1, ascov$Gamma, cov1)
fit1<-robust.fit1[[4]][[1]]
## model with the single factor
model<-paste('f =~ ', vname[1])
for (i in 2:length(vname)){
model<-paste(model, '+', vname[i])
}
## estimate the model using rsem
cfa.res2<-cfa(model, sample.cov=sigma, sample.mean=cov1$mu, meanstructure = TRUE, sample.nobs=n, std.lv=TRUE)
robust.fit2 <- rsem.fit(cfa.res2, ascov$Gamma, cov1)
fit2<-robust.fit2[[4]][[1]]
cat("Test of tau equivalent\n")
cat("The robust F statistic is ", round(fit1[1],3), "\n")
cat("with a p-value ", round(fit1[4],4), "\n\n")
cat("Test of homogeneous items\n")
cat("The robust F statistic is ", round(fit2[1],3), "\n")
cat("with a p-value ", round(fit2[4],4), "\n")
test<-list(tau=cfa.res1,homo=cfa.res2,tau.fit=fit1,homo.fit=fit2)
class(test)<-'test'
invisible(test)
}
bootstrap<-function(y, type="omega", alpha=.95, nboot=1000, ci="bc", plot=FALSE, varphi=0, complete=FALSE, auxiliary=NULL, silent=FALSE){
## save the results to a vector
boot.res<-rep(NA, nboot)
n<-nrow(y)
if (type=="omega"){
## for omega
temp<-try(omega(y,varphi=varphi, se=FALSE, test=FALSE, silent=silent))
orig.est<-temp$omega
for (i in 1:nboot){
bdata<-y[sample(1:n, n, replace=TRUE), ]
temp<-try(omega(bdata,varphi=varphi, se=FALSE, test=FALSE, silent=silent))
if (!inherits(temp, "try-error")) {
boot.res[i] <- temp$omega
}
}
}else{
## for alpha
temp<-try(alpha(y,varphi=varphi, se=FALSE, test=FALSE, silent=silent))
orig.est<-temp$alpha
for (i in 1:nboot){
bdata<-y[sample(1:n, n, replace=TRUE), ]
temp<-try(alpha(bdata,varphi=varphi, se=FALSE, test=FALSE, silent=silent))
if (!inherits(temp, "try-error")) {
boot.res[i] <- temp$alpha
}
}
}
boot.se<-sd(boot.res, na.rm=TRUE)
if (alpha>.5) alpha<-1-alpha
prob<-c(alpha/2, 1-alpha/2)
if (ci=="bc"){
z0<-qnorm(mean(boot.res<orig.est, na.rm=TRUE))
prob<-c(pnorm(2*z0+qnorm(alpha/2)), pnorm(2*z0+qnorm(1-alpha/2)))
boot.ci<-quantile(boot.res, prob=prob)
}else{
boot.ci<-quantile(boot.res, prob=prob)
}
## output the results
if (type=="omega"){
cat('The estimated omega is ', orig.est, "\n")
}else{
cat('The estimated alpha is ', orig.est, "\n")
}
cat('Its bootstrap se is ', round(boot.se, 3), "\n")
cat('Its bootstrap confidence interval is [', round(boot.ci[1], 3), " , ", round(boot.ci[2], 3),"]\n")
if (plot) plot(density(boot.res), xlab='Bootstrap distribution', ylab='density', main="")
bootobj<-list(estimate=orig.est, se=boot.se, ci=boot.ci, boot.est=boot.res, type=type)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.