qmboots <- function(dataset, nfactors, nsteps, load="auto", rotation="varimax", indet="qindtest", fsi=TRUE, forced=T, distribution=NULL, cor.method="pearson", ...) {
warning("Note that bootstrapping is an advanced technique in Q methodology, see the help page and relevant reference.")
if(forced == F) warning("Because the data are specified as non-forced, it's advisable to focus the interpretation on the z-scores (rather than on the factor scores)")
startime <- Sys.time()
nstat <- nrow(dataset)
nqsorts <- ncol(dataset)
# Number of iterations requested
itnumber <- paste("Number of requested iterations:", nsteps)
print(itnumber, quote=FALSE)
#-----------------------------------------------
# A. create objects and slots for the results
#-----------------------------------------------
#bootstrap results
qmbr <- vector("list", nfactors)
names(qmbr) <- paste0("factor_", 1:nfactors)
n <- 1
while (n <= nfactors) {
qmbr[[n]] <- list()
qmbr[[n]]$flagged <- data.frame(matrix(as.logical(NA), nrow=nqsorts, ncol=nsteps,
dimnames=list(colnames(dataset),
paste0("step_",1:nsteps))))
qmbr[[n]]$zsc <- data.frame(matrix(as.numeric(NA), nrow=nstat, ncol=nsteps,
dimnames=list(row.names(dataset),
paste0("step_",1:nsteps))))
qmbr[[n]]$loa <- data.frame(matrix(as.numeric(NA), nrow=nqsorts, ncol=nsteps,
dimnames=list(colnames(dataset),
paste0("step_",1:nsteps))))
n <- n+1
}
#indeterminacy tests results
if (indet == "procrustes") {
}
if (indet == "qindtest" | indet == "both") {
#create dataframe for results and reports from INDETERMINACY issue
qmts <- list()
qmts[[1]] <- data.frame(matrix(NA, nrow=nfactors, ncol=nsteps,
dimnames=list(paste0("f",c(1:nfactors)), paste0("order_",1:nsteps))))
qmts[[2]] <- data.frame(matrix(NA, nrow=nfactors, ncol=nsteps,
dimnames=list(paste0("f",c(1:nfactors)), paste0("sign_",1:nsteps))))
names(qmts) <- c("torder", "tsign")
qmts_log <- list()
qmts_log[[1]] <- data.frame(matrix(NA, nrow=1, ncol=nsteps, dimnames=list("log", paste0("log_ord_",1:nsteps))))
qmts_log[[2]] <- data.frame(matrix(NA, nrow=1, ncol=nsteps, dimnames=list("log", paste0("log_sig_",1:nsteps))))
names(qmts_log) <- c("torder", "tsign")
}
#-----------------------------------------------
# B. full sample Q method loadings (target matrix)
#-----------------------------------------------
#use manually introduced matrix of factor loadings as target, if any
if (is.matrix(load) | is.data.frame(load)) {
if ((nrow(load) == nqsorts) & (ncol(load) == nfactors)) {
flagged <- qflag(loa=load, nstat=nstat)
qm <- qzscores(dataset, nfactors, flagged=flagged, loa=load)
target <- load
colnames(target) <- paste0("target_f", 1:nfactors)
} else stop("Q method input: The factor loading matrix provided in 'load' should have the correct number of Q-sort as rows and the correct number of rotated factors as columns.")
} else if (is.character(load) & length(load) == 1) {
if (load == "auto") {
qm <- qmethod(dataset, nfactors, rotation=rotation, forced=forced, distribution=distribution, cor.method=cor.method, ...)
flagged <- qm$flagged
target <- as.matrix(qm$loa)
colnames(target) <- paste0("target_f", 1:nfactors)
} else stop("Q method input: 'load' has to be either 'auto' or a matrix.")
}
#create vector to build resamples
#list of number of Q sorts, repeated as many times as bootstrap replications, to ensure that each Q-sort appears the same amount of times
qsvector <- rep(1:nqsorts, times=nsteps)
#generate a random order for the previous list
qsrand <- sample(1:(nqsorts*nsteps), size=(nqsorts*nsteps), replace=FALSE)
#order qsvector randomly
qsvector <- qsvector[qsrand]
#-----------------------------------------------
# C. actual bootstrap
#-----------------------------------------------
it_count <- 1
while (it_count <= nsteps) {
qp <- 1 + (nqsorts*(it_count-1))
# Create bootstrap resample
subdata <- dataset[ , qsvector[c(qp:(qp+(nqsorts-1)))]]
# Reshape target matrix of original factor loadings
subtarget <- target[qsvector[c(qp:(qp+(nqsorts-1)))],]
# Full bootstrap step
step_res <- qbstep(subdata=subdata, subtarget=subtarget,
indet, nfactors, nqsorts, nstat,
qmts=qmts, qmts_log=qmts_log,
flagged=flagged, forced=forced,
distribution=distribution, cor.method=cor.method, ...)
# Export essential results: flagged, zsc and loa
for (n in 1:nfactors) {
# z-scores
qmbr[[n]][["zsc"]][paste0("step_",it_count)] <- step_res$zsc[[n]]
# Flagged Q-sorts and factor loadings -- important to assign to the correct columns!
for (r in 1:nqsorts) {
if (sum(rownames(qmbr[[n]][[3]])[r] == names(subdata)) != 0) { #if the Q sort is in the resample...
qmbr[[n]][["flagged"]][r,it_count] <- step_res$flagged[[n]][which(rownames(qmbr[[n]][["flagged"]])[r] == names(subdata))]
qmbr[[n]][["loa"]][r,it_count] <- step_res$loadings[[n]][which(rownames(qmbr[[n]][["loa"]])[r] == names(subdata))]
}
colnames(qmbr[[n]][["flagged"]])[it_count] <- paste0("step_",it_count)
colnames(qmbr[[n]][["loa"]])[it_count] <- paste0("step_",it_count)
}
}
# Export results of indeterminacy correction, if any
if (indet == "both" | indet == "qindtest") {
# Test results (logical)
qmts[[1]][paste("order_",it_count, sep="")] <- step_res[[4]]
qmts[[2]][paste("sign_",it_count, sep="")] <- step_res[[5]]
# Reports of solution implementation
qmts_log[[1]][paste("log_ord_",it_count, sep="")] <- step_res[[6]]
qmts_log[[2]][paste("log_sig_",it_count, sep="")] <- step_res[[7]]
}
# Count the iteration
it_msg <- paste("Finished iteration number", it_count)
print(it_msg, quote=FALSE)
it_count <- it_count + 1
}
#---actual bootstrap ends here
#-----------------------------------------------
# D. Export indeterminacy results into one object
#-----------------------------------------------
qindet <- list()
if (indet == "none") {
qindet <- "Caution: no correction of PCA bootstrap indeterminacy issue was performed. This may introduced inflated variability in the results"
}
if (indet == "procrustes") {
qindet <- "Procrustes rotation from MCMCpack was used to solve the PCA indeterminacy issue ('alignment problem')"
}
if (indet == "qindtest" | indet == "both") {
qindet[[1]] <- qmts
qindet[[2]] <- qmts_log
}
#-----------------------------------------------
# E. Test which steps could not be swap corrected
#-----------------------------------------------
if (indet == "qindtest" | indet == "both") {
errmsg <- "ERROR in ORDER swap: at least one factor in the resample is best match for two or more factors in the target"
badsteps <- which(qindet[[2]][[1]] == errmsg)
print(paste("Number of steps with order swap error:",
length(badsteps)))
if (sum(qindet[[2]][[1]] == errmsg) == 0) badsteps <- 0
} else badsteps <- 0
#-----------------------------------------------
# F. Summary stats for z-scores of bootstrap
#-----------------------------------------------
qmbs <- list() # Q method bootstrap summary
#- - - - - - - - - - - - - - - - - - - - - - - -
# 1. z-scores
for (n in 1:nfactors) {
if (sum(badsteps) > 0 & (indet == "qindtest" | indet == "both")) {
t.zsc <- qmbr[[n]]$zsc[,-badsteps]
} else t.zsc <- qmbr[[n]]$zsc
qmbs[[n+1]] <- merge(describe(t(t.zsc)),
t(apply(t.zsc, 1, quantile,
probs=c(0.025, 0.25, 0.75, 0.975),
na.rm=TRUE)),
by='row.names', sort=FALSE)
rownames(qmbs[[n+1]]) <- qmbs[[n+1]][,1]
qmbs[[n+1]][,1] <- NULL
}
names(qmbs) <- paste0("factor", 0:nfactors)
#- - - - - - - - - - - - - - - - - - - - - - - -
# 2. Factor scores (fragment adapted from qzscores.R)
if (forced==T) qscores <- sort(dataset[,1], decreasing=FALSE)
if (forced==F) qscores <- distribution
# Build frame for factor scores
zsc_mea <- data.frame(zsc_mea=c(1:nstat), row.names=row.names(dataset))
zsc_bn <- data.frame(zsc_bn=c(1:nstat), row.names=row.names(dataset))
for (f in 1:nfactors) zsc_mea[,f] <- qmbs[[f+1]]$mean
colnames(zsc_mea) <- paste0("zsc_mea_f", c(1:nfactors))
for (f in 1:nfactors) {
for (s in 1:nstat) {
# Find which statement has the current qscore rank
statement <- order(zsc_mea[,f])[[s]]
zsc_bn[statement,f] <- qscores[[s]]
}
}
colnames(zsc_bn) <- paste0("fsc_f", c(1:nfactors))
qmbs[[1]] <- zsc_bn
names(qmbs)[1] <- c("Bootstraped factor scores")
#-----------------------------------------------
# G. Summary stats for Q-sort factor loadings of bootstrap
#-----------------------------------------------
qmbl <- list()
for (n in 1:nfactors) {
if (sum(badsteps) > 0 & indet == "qindtest" | indet == "both") {
t.fla <- qmbr[[n]]$flagged[,-badsteps]
t.loa <- qmbr[[n]]$loa[,-badsteps]
} else {
t.fla <- qmbr[[n]]$flagged
t.loa <- qmbr[[n]]$loa
}
# Factor loadings
qmbl[[n]] <- merge(describe(t(t.loa)),
t(apply(t.loa, 1, quantile,
probs=c(0.025, 0.25, 0.75, 0.975),
na.rm=TRUE)),
by='row.names', sort=FALSE)
rownames(qmbl[[n]]) <- qmbl[[n]][,1]
qmbl[[n]][,1] <- NULL
# Frequency of flagging
qmbl[[n]]$flag_freq <- rowMeans(t.fla, na.rm=T)
for (j in 1:nqsorts) {
if (sum(badsteps) > 0 & indet == "qindtest" | indet == "both") {
token <- t(qmbr[[n]][[1]][,-badsteps])[,j]
} else token <- t(qmbr[[n]][[1]])[,j]
qmbl[[n]][j,"flag_freq"] <- length(which(token == TRUE)) / (length(which(token == TRUE)) + length(which(token == FALSE)))
}
qmbl[[n]] <- qmbl[[n]][order(row.names(qmbl[[n]])), ]
n <- n+1
}
names(qmbl) <- paste0("factor", 1:nfactors)
#-----------------------------------------------
# H. Export everything and report
#-----------------------------------------------
qmboots <- list()
qmboots[[1]] <- qmbs
qmboots[[2]] <- qmbr
qmboots[[3]] <- qindet
qmboots[[4]] <- as.data.frame(matrix(qsvector, nrow=nqsorts,
dimnames=list(NULL, paste0("bsampl_", 1:nsteps))))
qmboots[[5]] <- qm
qmboots[[6]] <- qscores
qmboots[[7]] <- qmbl
if (fsi == TRUE) {
# Calculate FACTOR STABILITY INDEX
fsii <- qfsi(nfactors=nfactors, nstat=nstat, qscores=qscores, zsc_bn=zsc_bn, qm=qm)
qmboots[[8]] <- fsii
names(qmboots) <- c("zscore-stats", "full.bts.res", "indet.tests", "resamples", "orig.res", "q.array", "loa.stats", "fsi")
} else {
names(qmboots) <- c("zscore-stats", "full.bts.res", "indet.tests", "resamples", "orig.res", "q.array", "loa.stats")
}
fintime <- Sys.time()
duration <- paste(format(floor(difftime(fintime, startime, units="hours")[[1]]), width=2),":", format(floor(difftime(fintime, startime, units="mins")[[1]] %% 60), width=2),":", format(floor(difftime(fintime, startime, units="secs")[[1]] %% 60), width=2), sep="")
cat("-----------------------------------------------\nBootstrap of ",nsteps," steps\nCall: qmboots(nfactors=",nfactors, ", nstat=", nstat, ", nqsorts=",nqsorts, ", nsteps=",nsteps, ", load=", load, ", rotation=", rotation, ", indet=",indet, ", fsi=",fsi,") \n-----------------------------------------------\nSTARTED : ",format(startime)," \nFINISHED: ", format(Sys.time()), "\nDURATION: ", duration, " hrs:min:sec","\n-----------------------------------------------\nAlignment correction method: ",indet,"\nNumber of steps that could not be reordered: ",length(badsteps), "\n", sep="")
invisible(qmboots)
}
#TO-DO: set printing method to show the bootstrapped z-scores, factor scores and stability indices
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.