Nothing
dimV <- function(x){
d <- dim(x)
if (is.null(d)){
d <- length(x)
}
d
}
# Wrapper for samples of fixed and samples of subjects:
parSamples_Mplus <- function(
fixed, # Array with fixed effect samples
SD = NULL, # Array with SD samples
subject = NULL
){
if (missing(fixed)){
stop("'fixed' must not be missing.")
}
# if (missing(SD)){
# stop("'SD' must not be missing.")
# }
aggDims <- 1:(length(dimV(fixed))-1)
mean <- apply(fixed, aggDims, mean)
lower <- apply(fixed, aggDims, quantile, 0.025)
upper <- apply(fixed, aggDims, quantile, 0.975)
if (!is.null(SD)){
SD <- apply(SD, aggDims, mean)
}
modelArray(
mean = mean,
lower = lower,
upper = upper,
SD = SD,
subject = subject
)
}
# Wrapper function to compute covs based on mplus samples
covSamples_Mplus <- function(fixed, SD = NULL, subject = NULL){
if (missing(fixed)){
stop("'fixed' must not be missing.")
}
# Covariance:
cov_mean <- apply(fixed, 1:2, mean)
cov_lower <- apply(fixed, 1:2, quantile, 0.025)
cov_upper <- apply(fixed, 1:2, quantile, 0.975)
cov_subject <- NULL
# cor samples:
corSamples <- fixed
for (i in 1:dimV(corSamples)[3]){
corSamples[,,i] <- cov2cor(corSamples[,,3])
}
# Correlation:
cor_mean <- apply(corSamples, 1:2, mean)
cor_lower <- apply(corSamples, 1:2, quantile, 0.025)
cor_upper <- apply(corSamples, 1:2, quantile, 0.975)
cor_subject <- NULL
# Precision samples:
precSamples <- fixed
for (i in 1:dimV(precSamples)[3]){
precSamples[,,i] <- corpcor::pseudoinverse(precSamples[,,i])
}
# Precision:
prec_mean <- apply(precSamples, 1:2, mean)
prec_lower <- apply(precSamples, 1:2, quantile, 0.025)
prec_upper <- apply(precSamples, 1:2, quantile, 0.975)
prec_subject <- NULL
# pcor:
pcorSamples <- fixed
for (i in 1:dimV(pcorSamples)[3]){
pcorSamples[,,i] <- corpcor::cor2pcor(pcorSamples[,,i])
}
pcor_mean <- apply(pcorSamples, 1:2, mean)
pcor_lower <- apply(pcorSamples, 1:2, quantile, 0.025)
pcor_upper <- apply(pcorSamples, 1:2, quantile, 0.975)
pcor_subject <- NULL
modelCov(
cov = modelArray(
mean = cov_mean,
lower = cov_lower,
upper = cov_upper,
subject = cov_subject
),
cor = modelArray(
mean = cor_mean,
lower = cor_lower,
upper = cor_upper,
subject = cor_subject
),
prec = modelArray(
mean = prec_mean,
lower = prec_lower,
upper = prec_upper,
subject = prec_subject
),
pcor = modelArray(
mean = pcor_mean,
lower = pcor_lower,
upper = pcor_upper,
subject = pcor_subject
)
)
}
# Wrapper for samples of fixed and samples of subjects:
parSamples_bugs <- function(
samples, # BUGS samples
fixed, # Name of fixed
subject, # Name of subject
SD # Name of SD samples
){
# Check if missing:
if (!missing(fixed)){
fixedSamples <- samples[[fixed]]
} else {
fixedSamples <- NULL
}
if (!missing(subject)){
subjectSamples <- samples[[subject]]
} else {
subjectSamples <- NULL
}
if (!missing(SD)){
SDSamples <- samples[[SD]]
} else {
SDSamples <- NULL
}
# Post-hoc estimate missing samples:
if (is.null(fixedSamples) && is.null(subjectSamples)){
stop("'fixed' and 'subject' are both missing")
}
# If fixed is missing, construct fixed from averaging subjects in each iteration:
if (is.null(fixedSamples) && !is.null(subjectSamples)){
dimSubject <- dimV(subjectSamples)
fixedSamples <- apply(subjectSamples,seq_len(length(dimSubject)-1),mean)
}
# If SD is missing, compute from subjects:
if (is.null(SDSamples) && !is.null(subjectSamples)){
dimSubject <- dimV(subjectSamples)
SDSamples <- apply(subjectSamples,seq_len(length(dimSubject)-1),sd)
}
# Subject list:
if (is.null(subjectSamples)){
subjectList <- NULL
} else {
dimSubject <- dimV(subjectSamples)
nDim <- length(dimSubject)
# Structure as list:
subjectMeans <- apply(subjectSamples, seq_len(nDim)[-1], mean)
subjectList <- list()
for (p in seq_len(dimSubject[nDim])){
if (nDim == 3){
# Vector
subjectList[[p]] <- subjectMeans[,p,drop=FALSE]
dim(subjectList[[p]]) <- dimSubject[2]
} else if (nDim == 4){
# Vector
subjectList[[p]] <- subjectMeans[,,p,drop=FALSE]
dim(subjectList[[p]]) <- dimSubject[2:3]
}else if (nDim == 5){
# Vector
subjectList[[p]] <- subjectMeans[,,,p,drop=FALSE]
dim(subjectList[[p]]) <- dimSubject[2:4]
} else stop(">3 dimensional arrays not supported")
}
}
aggDims <- 2:length(dimV(fixedSamples))
mean <- apply(fixedSamples, aggDims, mean)
lower <- apply(fixedSamples, aggDims, quantile, 0.025)
upper <- apply(fixedSamples, aggDims, quantile, 0.975)
if (!is.null(SDSamples)){
SD <- apply(SDSamples, aggDims, mean)
} else {
SD <- array(NA, dim = dimV(mean))
}
modelArray(
mean = mean,
lower = lower,
upper = upper,
SD = SD,
subject = subjectList
)
}
# Wrapper function to compute this based on samples of covariance:
covSamples_bugs <- function(samples, fixed, subject){
if (missing(subject)){
subjectSamples <- NULL
} else {
subjectSamples <- samples[[subject]]
}
# If fixed is missing, average from subject:
if (missing(fixed)){
if (missing(subject)){
stop("'subject' and 'fixed' may not both be missing.")
}
fixedSamples <- apply(subjectSamples,1:3,mean)
} else {
fixedSamples <- samples[[fixed]]
}
# Covariance:
cov_mean <- apply(fixedSamples, 2:3, mean)
cov_lower <- apply(fixedSamples, 2:3, quantile, 0.025)
cov_upper <- apply(fixedSamples, 2:3, quantile, 0.975)
if (!is.null(subjectSamples)){
cov_subject <- list()
nP <- dim(subjectSamples)[length( dim(subjectSamples))]
for (p in 1:nP){
cov_subject[[p]] <- apply(subjectSamples[,,,p],2:3,mean)
}
} else {
cov_subject <- NULL
}
# cor samples:
corSamples <- fixedSamples
for (i in 1:dimV(corSamples)[1]){
corSamples[i,,] <- cov2cor(corSamples[i,,])
}
# Correlation:
cor_mean <- apply(corSamples, 2:3, mean)
cor_lower <- apply(corSamples, 2:3, quantile, 0.025)
cor_upper <- apply(corSamples, 2:3, quantile, 0.975)
if (!is.null(subjectSamples)){
cor_subject <- list()
for (p in 1:nP){
corSamples <- subjectSamples[,,,p]
for (i in 1:dimV(corSamples)[1]){
corSamples[i,,] <- cov2cor(corSamples[i,,])
}
cor_subject[[p]] <- apply(corSamples,2:3,mean)
}
} else {
cor_subject <- NULL
}
# Precision samples:
precSamples <- fixedSamples
for (i in 1:dimV(precSamples)[1]){
precSamples[i,,] <- corpcor::pseudoinverse(precSamples[i,,])
}
# Precision:
prec_mean <- apply(precSamples, 2:3, mean)
prec_lower <- apply(precSamples, 2:3, quantile, 0.025)
prec_upper <- apply(precSamples, 2:3, quantile, 0.975)
if (!is.null(subjectSamples)){
prec_subject <- list()
for (p in 1:nP){
precSamples <- subjectSamples[,,,p]
for (i in 1:dimV(precSamples)[1]){
precSamples[i,,] <- corpcor::pseudoinverse(precSamples[i,,])
}
prec_subject[[p]] <- apply(precSamples,2:3,mean)
}
} else {
prec_subject <- NULL
}
# pcor:
pcorSamples <- fixedSamples
for (i in 1:dimV(pcorSamples)[1]){
pcorSamples[i,,] <- corpcor::cor2pcor(pcorSamples[i,,])
}
pcor_mean <- apply(pcorSamples, 2:3, mean)
pcor_lower <- apply(pcorSamples, 2:3, quantile, 0.025)
pcor_upper <- apply(pcorSamples, 2:3, quantile, 0.975)
if (!is.null(subjectSamples)){
pcor_subject <- list()
for (p in 1:nP){
pcorSamples <- subjectSamples[,,,p]
for (i in 1:dimV(pcorSamples)[1]){
pcorSamples[i,,] <- corpcor::cor2pcor(pcorSamples[i,,])
}
pcor_subject[[p]] <- apply(pcorSamples,2:3,mean)
}
} else {
pcor_subject <- NULL
}
modelCov(
cov = modelArray(
mean = cov_mean,
lower = cov_lower,
upper = cov_upper,
subject = cov_subject
),
cor = modelArray(
mean = cor_mean,
lower = cor_lower,
upper = cor_upper,
subject = cor_subject
),
prec = modelArray(
mean = prec_mean,
lower = prec_lower,
upper = prec_upper,
subject = prec_subject
),
pcor = modelArray(
mean = pcor_mean,
lower = pcor_lower,
upper = pcor_upper,
subject = pcor_subject
)
)
}
modelCov <- function(
cov,
cor, # Covariance
prec, # Precision (inverse cov)
pcor # Partial correlation
){
if (missing(cov)){
stop("'cov' can not be missing")
}
if (!missing(cor) && !is(cor,"mlVARarray")){
stop("'cor' must be missing or an object of class 'mlVARarray'")
}
if (!missing(prec) && !is(prec,"mlVARarray")){
stop("'prec' must be missing or an object of class 'mlVARarray'")
}
if (!missing(pcor) && !is(pcor,"mlVARarray")){
stop("'pcor' must be missing or an object of class 'mlVARarray'")
}
cov2corNA <- function(x){
if (any(is.na(x))){
x[] <- NA
return(x)
} else {
return(cov2cor(x))
}
}
solveNA <- function(x){
if (any(is.na(x))){
x[] <- NA
return(x)
} else {
return(corpcor::pseudoinverse(x))
}
}
cor2pcorNA <- function(x){
if (any(is.na(x))){
x[] <- NA
return(x)
} else {
return(cor2pcor(x))
}
}
# Fill missings (only mean and subject):
if (missing(cor)){
cor <- modelArray(
mean = cov2corNA(cov[['mean']]),
subject = lapply(cov[['subject']],cov2corNA)
)
}
if (missing(prec)){
prec <- modelArray(
mean = solveNA(cov[['mean']]),
subject = lapply(cov[['subject']],solveNA)
)
}
if (missing(pcor)){
pcor <- modelArray(
mean = cor2pcorNA(cov[['mean']]),
subject = lapply(cov[['subject']],cor2pcorNA)
)
}
Results <- list(
cov=cov,
cor=cor,
prec=prec,
pcor=pcor
)
class(Results) <- c("mlVarCov","list")
return(Results)
}
modelArray <- function(
mean, # If missing computed from subject
subject, # LIST containing estimate per subject
SE, # Can be missing
P, # If missing computed from mean and SE
lower, # If missing computed from mean and SE
upper,
SD
){
Results <- list()
if (missing(subject)){
subject <- NULL
}
# Either mean or subject must not be missing:
if (is.null(subject) && missing(mean)){
stop("Either 'subject' or 'mean' must be used")
}
# If mean missing but not subject, compute:
if (missing(mean) && !is.null(subject)){
N <- length(subject)
mean <- Reduce("+",subject) / N
}
# Dimensions:
dim <- dimV(mean)
if (is.null(dim)){
dim <- length(mean)
}
# If SD missing but not subject, compute:
if (missing(SD)){
if (!is.null(subject)){
N <- length(subject)
mean_for_sd <- Reduce("+",subject) / N
SD <- sqrt(Reduce('+',lapply(subject,function(x)(x - mean_for_sd)^2)) / (N-1))
} else {
SD <- array(NA, dim)
}
}
# If P is missing but SE is not, compute:
if (missing(P)){
if (!missing(SE)){
P <- 2*(1-pnorm(abs(mean/SE)))
} else {
P <- array(NA, dim)
}
}
# If SE is missing, NA:
if (missing(SE)){
SE <- array(NA, dim)
}
# If lower is missing but mean and SE not, compute:
if (missing(lower)){
if (!missing(SE)){
lower <- mean - 1.959964 * SE
} else {
lower <- array(NA, dim)
}
}
# If upper is missing but mean and SE not, compute:
if (missing(upper)){
if (!missing(SE)){
upper <- mean + 1.959964 * SE
} else {
upper <- array(NA, dim)
}
}
# Subject:
if (is.null(subject)){
subject <- NULL
}
# Store mean to list:
Results[['mean']] <- mean
Results[['SD']] <- SD
Results[['lower']] <- lower
Results[['upper']] <- upper
Results[['SE']] <- SE
Results[['P']] <- P
Results[['subject']] <- subject
class(Results) <- c("mlVARarray","list")
return(Results)
}
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.