Nothing
# helpers.R
# compute interaction, squared, cubic, dummy variables, etc. for RSA
add.variables <- function(formula, df) {
IV1 <- all.vars(formula)[2]
IV2 <- all.vars(formula)[3]
IV12 <- paste0(IV1, "2")
IV22 <- paste0(IV2, "2")
IV13 <- paste0(IV1, "3")
IV23 <- paste0(IV2, "3")
IV_IA <- paste0(IV1, "_", IV2)
IV_IA2 <- paste0(IV1, "2", "_", IV2)
IV_IA3 <- paste0(IV1, "_", IV2, "2")
df[, IV12] <- df[, IV1]^2
df[, IV22] <- df[, IV2]^2
df[, IV_IA] <- df[, IV1]*df[, IV2]
# three new variables for piecewise regression (test absolute difference score) - Edwards (2002) model
df$W.JRE <- ifelse(df[, IV1] >= df[, IV2], 0, 1)
df[, paste0("W.JRE_", IV1)] <- df$W.JRE*df[, IV1]
df[, paste0("W.JRE_", IV2)] <- df$W.JRE*df[, IV2]
# three new variables for piecewise regression (test absolute difference score) - new model Schoenbrodt 2012
df$W <- ifelse(df[, IV1] >= df[, IV2], 1, -1)
df$W[df[, IV1] == df[, IV2]] <- 0
df[, paste0("W_", IV1)] <- df$W*df[, IV1]
df[, paste0("W_", IV2)] <- df$W*df[, IV2]
df$diff <- df[, IV2] - df[, IV1]
df$SD <- df$diff^2
df$absdiff <- abs(df$diff)
# cubic terms
df[, IV13] <- df[, IV1]^3
df[, IV_IA2] <- df[, IV1]^2*df[, IV2]
df[, IV_IA3] <- df[, IV1]*df[, IV2]^2
df[, IV23] <- df[, IV2]^3
return(df)
}
# helper function: takes a list of lavaan models (can include NULLs), and returns the usual anova object
anovaList <- function(modellist) {
mods <- modellist[!sapply(modellist, function(x) is.null(x))]
mods <- mods[!sapply(mods, function(x) !inspect(x, "converged"))]
if (length(mods) == 0) {
return(list(n.mods=0))
}
# put them in order (using df)
DF <- sapply(mods, fitmeasures, "df")
mods <- mods[order(DF, decreasing = FALSE)]
# prevent lavaan error in case that some models (but not all) have scaled test statistics
# detect such cases (condition copy-pasted from lav_test_LRT.R in lavaan)
mods.scaled <- unlist( lapply(mods, function(x) {
any(c("satorra.bentler", "yuan.bentler", "yuan.bentler.mplus", "mean.var.adjusted", "scaled.shifted")
%in% unlist(sapply(slot(x, "test"), "[", "test")) ) }))
# change internal label
if( ! ( all(mods.scaled) | !any(mods.scaled) ) ) {
mods[[ which(sapply(mods, fitmeasures, "df") == 0) ]]@test[[2]]$test <- mods[[ which(mods.scaled)[1] ]]@test[[2]]$test
}
pStr <- sapply(1:length(mods), function(x){
if(x==1) {
paste("mods[[",x,"]]",sep = "")
} else {
paste("force(mods[[",x,"]])",sep = "")
}
})
#pStr2 <- paste0("lavTestLRT(", paste(pStr, collapse=", "), ", method='satorra.bentler.2010')")
pStr2 <- paste0("lavTestLRT(", paste(pStr, collapse=", "), ", method='default')")
a1 <- eval(parse(text = pStr2))
if (length(mods) > 1) {
rownames(a1) <- names(mods)
}
attr(a1, "n.mods") <- length(mods)
return(list(ANOVA=a1, models=mods, n.mods=length(mods)))
}
## internal helper function: compare models
# mL = model list
# set = label that is attached to the results
cModels <- function(mL, set, free.max) {
aL1 <- anovaList(mL)
if (aL1$n.mods > 1) {
N <- lavaan::nobs(aL1$models[[1]])
a1 <- cbind(aL1$ANOVA[, c(1, 4:7)], plyr::ldply(aL1$models, function(X) {
F <- fitmeasures(X)
R <- inspect(X, "r2")
names(R) <- "R2"
n <- lavaan::nobs(X)
k <- free.max - F["df"] # number of parameters, including coefficients of control variables (if there are any)
suppressWarnings({
R2.p <- ifelse(k==0,
NA,
pf(((n-k-1)*R)/(k*(1-R)), k, n-k-1, lower.tail=FALSE))
})
names(R2.p) <- "R2.p"
# compute AICc
K <- k + 2 # number of parameters, including coefficients of control variables (if there are any), intercept and residual variance (thus the +2)
AICc <- -2*F["logl"] + 2*K + 2*(K*(K+1))/(n-K-1)
names(AICc) <- NULL
return(c(AICc=AICc, F[c("cfi", "srmr")], R, R2.p))
}))
a1 <- a1[, !grepl(".id", colnames(a1))]
a1$k <- free.max - a1$Df
a1$R2.adj <- 1 - ((1-a1$R2))*((N-1)/(N-a1$k-1))
a1$delta.R2 <- c(NA, a1$R2[1:(nrow(a1)-1)] - a1$R2[2:(nrow(a1))])
a1$model <- rownames(a1)
a1$set <- set
return(a1)
}
}
# internal helper function: F-test to compare R^2 difference between two nested models
# x = output object of RSA()
# unrestricted = name of unrestricted model (model with more estimated parameters)
# restricted = name of restricted model (less estimated parameters); with default setting "interceptonly", the function will return R^2 and respective pvalue of the unrestricted model
R2difftest <- function(x, unrestricted="", restricted="interceptonly"){
n <- lavaan::nobs(x$models[[unrestricted]])
free.max <- getFreeParameters(x$models[[unrestricted]])
Fu <- fitmeasures(x$models[[unrestricted]])
Ru <- inspect(x$models[[unrestricted]], "r2")
ku <- free.max - Fu["df"]
if ( restricted=="interceptonly" ){
Rr <- 0
kr <- 0
} else {
Fr <- fitmeasures(x$models[[restricted]])
Rr <- inspect(x$models[[restricted]], "r2")
kr <- free.max - Fr["df"]
}
suppressWarnings({
R2.p <- ifelse(ku==kr,
NA,
pf( ( ( n - ku - 1 ) * ( Ru - Rr ) ) / ( (ku-kr) * (1 - Ru ) ), ku-kr, n-ku-1, lower.tail=FALSE)
)
})
delta.R2 <- Ru-Rr
names(delta.R2) <- "R2"
names(R2.p) <- "R2.p"
return(list(
delta.R2=delta.R2,
R2.p=R2.p,
ku=ku
))
}
# simple wrapper: formats a number in f.2 format
f2 <- function(x, digits=2, prepoint=0, skipZero=FALSE) {
if (skipZero == TRUE) {zero <- "."} else {zero <- "0."}
if (length(dim(x)) == 2) {
apply(x, 2, function(x2) {gsub("0.", zero, sprintf(paste("%",prepoint,".",digits,"f",sep=""), x2) , fixed=TRUE)})
} else {
gsub("0.", zero, sprintf(paste("%",prepoint,".",digits,"f",sep=""), x) , fixed=TRUE)
}
}
# converts p values in stars
p2star <- function(val) {
res <- val
for (i in 1:length(val)) {
res[i] <- ""
if (is.na(val[i])) next();
if (val[i] <= 0.1) res[i] <- "\U2020"
if (val[i] <= 0.05) res[i] <- "*"
if (val[i] <= 0.01) res[i] <- "**"
if (val[i] <= 0.001) res[i] <- "***"
}
return(res)
}
# nicely formats a p-value
p0 <- function(x) {
if (is.na(x)) return("NA")
if (x >= .001) return(paste0("p = ", f2(x, 3, skipZero=TRUE)))
if (x < .001) return("p <.001")
}
p <- Vectorize(p0)
# returns number of maximum free parameters of a regression model
getFreeParameters <- function(model) {
VARS <- nrow(inspect(model, "free")$beta) # number of variables
df.max <- (VARS*(VARS+1))/2 # maximum df
df.pred <- ((VARS-1)*(VARS))/2 + 1 # df bound in the predictors (i.e., (co)variances of the predictors & variance of DV)
free.max <- df.max - df.pred # maximum of free parameters
return(free.max)
}
# computes the coordinates of an arbitrary intersection of the surface,
# defined by a line on the X-Y plane (p0 = intercept, p1=slope)
getIntersect <- function(b0=0, x=0, y=0, x2=0, xy=0, y2=0, p0, p1, xlim=c(-2, 2), grid=21) {
X <- seq(min(xlim), max(xlim), length.out=grid)
Y <- p0 + p1*X
n <- data.frame(X, Y)
n2 <- add.variables(z~X+Y, n)
n2$Z <- b0 + colSums(c(x, y, x2, y2, xy)*t(n2[, c(1:5)]))
return(n2[, c("X", "Y", "Z")])
}
model <- function(x, model="full") x$models[[model]]
syntax <- function(x, model="full") cat(x$models[[model]]@Options$model)
# transforms p-values to colors
pRamp <- function(p, sig=.05, borderline=.10, bias=.8) {
# calculate bias that the color transition is at the borderline value
bias2 <- .33/(borderline/(1 - sig))
cR1 <- colorRamp(c("red", "red", "orange"), bias=bias, space="Lab")
cR2 <- colorRamp(c("orange", "green", "green"), bias=bias2, space="Lab")
p2 <- rep("#FFFFFF", length(p))
if (length(p[p < sig])>0) {
p2[p < sig] <- rgb(cR1(p[p < sig]/sig), maxColorValue=255)
}
if (length(p[p >= sig])>0) {
p2[p >= sig] <- rgb(cR2((p[p >= sig] - sig) / (1 - sig)), maxColorValue=255)
}
return(p2)
}
# helper function: find closest value in vector
f0 <- function (vec, target, unique = TRUE) {
ret <- vec[sapply(target, function(x) which.min(abs(x - vec)))]
if (unique) { ret <- unique(ret) }
ret
}
# compute the predicted value from a single pair of predictors
predictRSA <- function(object, X, Y, model="full") {
C <- coef(object$models[[model]])
if (object$models[[model]]@Options$estimator != "DWLS") {
b0 <- as.numeric(ifelse(is.na(C[paste0(object$DV, "~1")]), b0, C[paste0(object$DV, "~1")]))
} else {
# the threshold is the negative of the intercept ...
b0 <- -as.numeric(ifelse(is.na(C[paste0(object$DV, "|t1")]), b0, C[paste0(object$DV, "|t1")]))
}
x <- as.numeric(ifelse(is.na(C["b1"]), 0, C["b1"]))
y <- as.numeric(ifelse(is.na(C["b2"]), 0, C["b2"]))
x2 <- as.numeric(ifelse(is.na(C["b3"]), 0, C["b3"]))
y2 <- as.numeric(ifelse(is.na(C["b5"]), 0, C["b5"]))
xy <- as.numeric(ifelse(is.na(C["b4"]), 0, C["b4"]))
w <- as.numeric(ifelse(is.na(C["w1"]), 0, C["w1"]))
wx <- as.numeric(ifelse(is.na(C["w2"]), 0, C["w2"]))
wy <- as.numeric(ifelse(is.na(C["w3"]), 0, C["w3"]))
# cubic parameters
x3 <- as.numeric(ifelse(is.na(C["b6"]), 0, C["b6"]))
x2y <- as.numeric(ifelse(is.na(C["b7"]), 0, C["b7"]))
xy2 <- as.numeric(ifelse(is.na(C["b8"]), 0, C["b8"]))
y3 <- as.numeric(ifelse(is.na(C["b9"]), 0, C["b9"]))
C <- c(x, y, x2, y2, xy, w, wx, wy,x3, x2y, xy2, y3)
# compute predicted value
Z <- b0 + colSums(C*t(cbind(X, Y, X^2, Y^2, X*Y, 0, 0, 0, X^3, X^2*Y, X*Y^2, Y^3)))
return(Z)
}
# fills up the long edges of a polygon with intermediate points
# If an edge is longer than minDist, new points re inserted.
# @param x Vector of x values
# @param y Vector of y values
interpolatePolygon <- function(x, y, minDist, plot=FALSE) {
minDist <- minDist^2 # compare with squared x^2 + y^2 (faster)
interp <- data.frame()
pol <- data.frame(x, y)
colnames(pol) <- c("x", "y")
for (i in 1:(nrow(pol)-1)) {
# get distance
D <- (pol[i, 1] - pol[i+1, 1])^2 + (pol[i, 2] - pol[i+1, 2])^2
if (D > minDist) {
N <- ceiling(sqrt(D)/sqrt(minDist)) # number of interpolations
APPROX <- data.frame(
x = seq(pol[i, 1], pol[i+1, 1], length.out=N),
y = seq(pol[i, 2], pol[i+1, 2], length.out=N)
)
interp <- rbind(interp, APPROX)
} else if (D>0 & D <= minDist){
interp <- rbind(interp, pol[i, ])
if (i==1) colnames(interp) <- c("x", "y")
}
}
interp <- rbind(interp, pol[nrow(pol), ])
if (plot==TRUE) {
plot(pol, col="red")
points(interp[, 1], interp[, 2], col="green", pch=20)
lines(interp[, 1], interp[, 2], col="darkgreen")
text(pol, label=1:nrow(pol))
}
return(interp)
}
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.