The MBESS R package, now on GitHub beginning with version 4.8.0.

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

DESCRIPTION

@@ -0,0 +1,26 @@
+Package: MBESS
+Type: Package
+Title: The MBESS R Package
+Version: 4.8.0
+Date: 2020-08-04
+Authors@R: c(person("Ken", "Kelley", role=c("aut", "cre"), email="kkelley@nd.edu"))
+Maintainer: Ken Kelley <kkelley@nd.edu>
+Depends: R (>= 3.5.0), stats
+Imports: boot, gsl, lavaan, MASS, methods, mnormt, nlme, OpenMx, parallel, sem, semTools 
+Description: Implements methods that useful in designing research studies and analyzing data, with 
+	particular emphasis on methods that are developed for or used within the behavioral, 
+	educational, and social sciences (broadly defined). That being said, many of the methods 
+	implemented within MBESS are applicable to a wide variety of disciplines. MBESS has a 
+	suite of functions for a variety of related topics, such as effect sizes, confidence intervals 
+	for effect sizes (including standardized effect sizes and noncentral effect sizes), sample size
+	planning (from the accuracy in parameter estimation [AIPE], power analytic, equivalence, and 
+	minimum-risk point estimation perspectives), mediation analysis, various properties of 
+	distributions, and a variety of utility functions. MBESS (pronounced 'em-bes') was originally 
+	an acronym for 'Methods for the Behavioral, Educational, and Social Sciences,' but at this 
+	point MBESS contains methods applicable and used in a wide variety of fields and is an 
+	orphan acronym, in the sense that what was an acronym is now literally its name. MBESS has 
+	greatly benefited from others, see <http://nd.edu/~kkelley/site/MBESS.html> for a detailed 
+	list of those that have contributed and other details.
+License: GPL-2 | GPL-3
+URL: http://nd.edu/~kkelley/site/MBESS.html
+RoxygenNote: 6.0.0

MBESS.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

NAMESPACE

@@ -0,0 +1,13 @@
+exportPattern("^[^\\.]") # Exports all functions not beginning with a period.
+import(methods)
+importFrom("grDevices", "dev.off", "jpeg", "pdf", "postscript",
+             "trans3d")
+importFrom("graphics", "abline", "arrows", "legend", "lines", "par",
+             "persp", "plot", "points", "rect", "segments", "text",
+             "title")
+importFrom("stats", "anova", "coef", "cor", "cov", "cov2cor", "fitted",
+             "integrate", "lm", "median", "na.omit", "nlm", "optim",
+             "optimize", "pchisq", "pf", "pnorm", "pt", "qchisq", "qf",
+             "qnorm", "qt", "quantile", "resid", "rnorm", "sd", "var",
+             "weighted.mean", "model.frame")
+importFrom("utils", "capture.output", "read.table", "write.table")

R/CFA.1.R

@@ -0,0 +1,108 @@
+CFA.1 <- function(S, N, equal.loading = FALSE, equal.error = FALSE, package = "lavaan", se = "standard", ...) 
+{
+if (!isSymmetric(S, tol = 1e-05)) stop("Input a symmetric covariance or correlation matrix 'S'")
+    q <- nrow(S)
+
+
+    if (package == "sem") 
+    {
+      if(!requireNamespace("sem", quietly = TRUE)) stop("The package 'sem' is needed; please install the package and try again (or use set 'package' to 'lavaan'.")
+
+        x <- matrix(NA, nrow = q, ncol = 1)
+        x <- paste("x", row(x), sep = "")
+
+        if (equal.loading) {
+            lamda <- matrix(rep("lamda", q), nrow = q, ncol = 1)
+        } else {
+            lamda <- matrix(NA, nrow = q, ncol = 1)
+            lamda <- paste("lamda", row(lamda), sep = "")
+        }
+
+        if (equal.error) {
+            psi.sq <- matrix(rep("psi.sq", q), nrow = q, ncol = 1)
+        } else {
+            psi.sq <- matrix(NA, nrow = q, ncol = 1)
+            psi.sq <- paste("psi.sq", row(psi.sq), sep = "")
+        }
+
+        model.1 <- cbind(paste("ksi", "->", x), lamda)
+        model.2 <- cbind(paste(x, "<->", x), psi.sq)
+        model <- rbind(model.1, model.2)
+        start <- matrix(NA, nrow = nrow(model), ncol = 1)
+        model <- cbind(model, start)
+        model <- rbind(model, c(paste("ksi", "<->", "ksi"), NA, 1))
+        class(model) <- "mod"
+
+        rownames(S) <- x
+        colnames(S) <- x
+        model.fitted <- sem::sem(model, S, N)
+        converged <- model.fitted$convergence
+        if (converged) {
+            if (equal.loading) 
+                k <- 1 else k <- q
+            Factor.Loadings <- model.fitted$coeff[1:k]
+            if (equal.error) 
+                m <- 1 else m <- q
+            Indicator.var <- model.fitted$coeff[k + 1:m]
+            Parameter.cov <- model.fitted$vcov
+        } else {
+            Factor.Loadings <- NA
+            Indicator.var <- NA
+            Parameter.cov <- NULL
+        }
+        result <- list(Model = model, Factor.Loadings = Factor.Loadings, Indicator.var = Indicator.var, Parameter.cov = Parameter.cov, converged = converged, package="sem")
+        return(result)
+    } 
+
+if (package == "lavaan") 
+  {
+if(!requireNamespace("lavaan", quietly = TRUE)) stop("The package 'lavaan' is needed; please install the package and try again.")
+
+        colnames(S) <- rownames(S) <- paste("y", 1:q, sep = "")
+
+        if (equal.loading) {
+            loadingName <- rep("a1", q)
+        } else {
+            loadingName <- paste("a", 1:q, sep = "")
+        }
+
+        if (equal.error) {
+            errorName <- rep("b1", q)
+        } else {
+            errorName <- paste("b", 1:q, sep = "")
+        }
+        model <- "f1 =~ NA*y1 + "
+        loadingLine <- paste(paste(loadingName, "*", colnames(S), sep = ""), collapse = " + ")
+        factorLine <- "f1 ~~ 1*f1\n"
+        errorLine <- paste(paste(colnames(S), " ~~ ", errorName, "*", colnames(S), 
+            sep = ""), collapse = "\n")
+        model <- paste(model, loadingLine, "\n", factorLine, errorLine, "\n")
+        try(fit <- lavaan::lavaan(model, sample.cov = S, sample.nobs = N, se = se, ...), silent = TRUE)
+
+        converged <- fit@Fit@converged
+        loading <- unique(as.vector(fit@Model@GLIST$lambda))
+        error <- unique(diag(fit@Model@GLIST$theta))
+        if (!all(error > 0)) 
+            converged <- FALSE
+        if (converged) {
+            if (se == "none") {
+                paramCov <- NULL
+            } else {
+                paramCov <- lavaan::vcov(fit)
+            }
+        } else {
+            loading <- NA
+            error <- NA
+            if (se == "none") {
+                paramCov <- NULL
+            } else {
+                paramCov <- NA
+            }
+        }
+        print("here")
+        result <- list(Model = model, Factor.Loadings = loading, Indicator.var = error, 
+            Parameter.cov = paramCov, converged = converged, package="lavaan")
+
+} 
+    return(result)
+} 

R/Expected.R2.R

@@ -0,0 +1,9 @@
+Expected.R2 <- function(Population.R2, N, p)
+{
+if(!requireNamespace("gsl", quietly = TRUE)) stop("The package 'gsl' is needed; please install the package and try again.")  
+
+Value <- 1 - ((N-p-1)/(N-1))*(1-Population.R2)*gsl::hyperg_2F1(1, 1, .5*(N+1), Population.R2)
+Value <- max(0, Value)
+Value <- min(Value, 1)
+return(Value)
+}

R/F2Rsquare.R

@@ -0,0 +1,6 @@
+"F2Rsquare" <-
+function(F.value=NULL, df.1=NULL, df.2=NULL)
+{
+if(is.null(df.1) | is.null(df.2)) stop("You have not specified \'df.1\' and/or \'df.2\'.")
+return(F.value*df.1/(F.value*df.1 + df.2))
+}

R/Lambda2Rsquare.R

@@ -0,0 +1,6 @@
+"Lambda2Rsquare" <-
+function(Lambda=NULL, N=NULL)
+{
+if(is.null(Lambda) | is.null(N)) stop("You must specify \'Lambda\' (i.e., a noncentral parameter) and \'N\' (i.e., sample size) in order to calculate the noncentraility paramter.")
+return(Lambda/(Lambda + N))
+}

R/Rsquare2F.R

@@ -0,0 +1,11 @@
+"Rsquare2F" <-
+function(R2=NULL, df.1=NULL, df.2=NULL, p=NULL, N=NULL)
+{
+if(is.null(df.1) & is.null(df.2) & !is.null(N) & !is.null(p))
+{
+df.1 <- p
+df.2 <- N-p-1
+}
+if(is.null(df.1) | is.null(df.2)) stop("You have not specified \'df.1\', \'df.2\', \'N\', and/or \'p\' correctly.")
+return((R2/df.1)/((1-R2)/df.2))
+}

R/Rsquare2Lambda.R

@@ -0,0 +1,6 @@
+"Rsquare2Lambda" <-
+function(R2=NULL, N=NULL)
+{
+if(is.null(R2) | is.null(N)) stop("You must specify \'R2\' (i.e., R Square) and \'N\' (i.e., sample size) in order to calculate the noncentraility paramter.")
+return((R2/(1-R2))*N)
+}

R/Sigma.2.SigmaStar.R

@@ -0,0 +1,145 @@
+Sigma.2.SigmaStar<- function(model, model.par, latent.var, discrep, ML=TRUE){
+
+if(is.null(names(model.par))) stop("The elements in 'model.par' must have names")
+if(sum(is.na(names(model.par))) !=0) stop ("Some of the elements in 'model.par' do not have names")
+if(length(unique(names(model.par))) != length(model.par)) stop("More than one element in 'model.par' has the same name")
+
+duplication.matrix <- function (n = 1) 
+{
+    if ((n < 1) | (round(n) != n)) 
+        stop("n must be a positive integer")
+    d <- matrix(0, n * n, n * (n + 1)/2)
+    count = 0
+    for (j in 1:n) {
+        d[(j - 1) * n + j, count + j] = 1
+        if (j < n) {
+            for (i in (j + 1):n) {
+                d[(j - 1) * n + i, count + i] <- 1
+                d[(i - 1) * n + j, count + i] <- 1
+            }
+        }
+        count = count + n - j
+    }
+    return(d)
+}
+
+vech<-function (x) 
+{
+    if (!is.square.matrix(x)) 
+        stop("argument x is not a square numeric matrix")
+    return(t(t(x[!upper.tri(x)])))
+}
+
+is.square.matrix<-function (x) 
+{
+    if (!is.matrix(x)) 
+        stop("argument x is not a matrix")
+    return(nrow(x) == ncol(x))
+}
+
+matrix.trace<-function (x) 
+{
+    if (!is.square.matrix(x)) 
+        stop("argument x is not a square matrix")
+    return(sum(diag(x)))
+}
+
+t<- length(model.par)
+r<- length(latent.var)
+T <- rep(0, t)
+R <- rep(0, r)
+#for(i in 1:t){
+#    T[i]<- sum(na.omit(unique(model[,2])==names(model.par)[i]))
+#    }
+#if(sum(T) != t) stop ("Some elements in 'model.par' have not been included in 'model'")
+
+delta<-discrep
+gamma0<- model.par
+res<-theta.2.Sigma.theta(model=model, theta=gamma0, latent.vars=latent.var)
+Sigma.gamma0<- res$Sigma.theta
+q <- res$t  # No. of model parameters
+p <- res$n  # No. of manifest variables
+p.star <- p*(p+1)/2
+
+D.mat <- duplication.matrix(p)
+#D.mat.MPinv <- invgen(D.mat)
+#K <- t(D.mat.MPinv)
+D <- t(D.mat)%*%D.mat
+if(isTRUE(ML)) W <- Sigma.gamma0
+W.inv <- solve(W)
+
+h<- 1e-8
+Sigma.deriv<- array(NA, c(p,p,q))
+B <- matrix(NA, p.star, q)
+for (i in 1:q){
+	u <- matrix(0, q,1)
+	u[i,]<-1
+	gamma<- gamma0+u*h
+	names(gamma)<-names(gamma0)
+	res.h<-theta.2.Sigma.theta(model=model, theta=gamma, latent.vars=latent.var)
+	Sigma.gamma <- res.h$Sigma.theta
+	Sigma.deriv[,,i] <- 	(Sigma.gamma-Sigma.gamma0)*(1/h)
+	B[,i]<- (-1)*D%*%vech(W.inv%*%Sigma.deriv[,,i]%*%W.inv) 
+	}
+
+y <- matrix(1:p.star/100, p.star,1)
+B.qr<- qr(B)
+e.tilt <- qr.resid(B.qr, y)
+
+E1 <- matrix(0, p,p)
+index<-1
+for (i2 in 1:p){
+	for(i1 in i2:p){
+		E1[i1, i2]<- e.tilt[index,1]
+		index<-index+1
+		}
+	}
+
+E2 <- matrix(0, p,p)
+index <- 1
+for (i1 in 1:p){
+  for (i2 in i1:p){
+    E2[i1, i2]<- e.tilt[index, 1]
+    index <- index+1
+    }
+  }
+
+E.tilt <- E1+E2-diag(diag(E1))	
+
+#E.tilt[upper.tri(E.tilt)]<- E.tilt[lower.tri(E.tilt)]
+#E.tilt[upper.tri(E.tilt)]<-0
+
+#E.tilt <- matrix(0, p,p)
+#E.tilt[lower.tri(E.tilt, diag=TRUE)]<- e.tilt
+#E.diag<- diag(E.tilt)
+#E.tilt <- E.tilt+t(E.tilt)-E.diag
+
+G <- W.inv %*% E.tilt
+get.kappa <- function(kappa, G, I, delta){
+	target<-abs(kappa*matrix.trace(G) - log(det(I+kappa*G))-delta)
+	return(target)
+	}
+
+kappa0 <- sqrt(2*delta/matrix.trace(G%*%G))
+I <- diag(p)
+res.kappa<- suppressWarnings(nlm(get.kappa, kappa0, G=G, I=I, delta=delta))
+kappa <- res.kappa$estimate
+iter<- res.kappa$iterations
+
+kappa<-as.numeric(kappa)
+E <- kappa*E.tilt
+Sigma.star <-Sigma.gamma0+E
+
+result <- list()
+result$Sigma.star<- Sigma.star
+result$Sigma_theta <-Sigma.gamma0
+#result$B <-B
+#result$E.tilt <-E.tilt
+result$E <-E
+#result$kappa0 <- kappa0
+#result$kappa <- kappa
+#result$iter <- iter
+#result$diag.E1<- diag(E1)
+#result$diag.E2 <- diag(E2)
+return(result)
+}#end of Sigma.2.SigmaStar<- function()

R/Variance.R2.R

@@ -0,0 +1,6 @@
+Variance.R2 <- function(Population.R2, N, p)
+{
+if(!requireNamespace("gsl", quietly = TRUE)) stop("The package 'gsl' is needed; please install the package and try again.")  
+
+return((((N-p-1)*(N-p+1))/(N^2-1))*((1-Population.R2)^2)*(gsl::hyperg_2F1(2, 2, .5*(N+3), Population.R2)) - ((Expected.R2(Population.R2=Population.R2, N=N, p=p) - 1)^2))
+}