Nothing
varPart <-
function(A, B, C = NA, AB, AC = NA, BC = NA, ABC = NA, model.type = NULL,
A.name = "Factor A", B.name = "Factor B", C.name = "Factor C",
model = NULL, groups = NULL, pred.type = "Y", cor.method = "pearson",
return.models = FALSE, plot = TRUE, plot.digits = 3, cex.names = 1.5,
cex.values = 1.2, main = "", cex.main = 2, plot.unexpl = TRUE,
colr = FALSE) {
# version 2.1 (21 Jul 2022)
if (!is.null(model.type)) message ("NOTE: Argument 'model.type' is no longer used.")
if (!is.null(model)) {
if (!missing(A) || !missing(B)) warning("Arguments 'A', 'B', etc. are ignored as argument 'model' is provided.")
if (!(inherits(model, "glm"))) stop("Argument 'model' is currently only implemented for class 'glm'.")
if (!(pred.type %in% c("Y", "P", "F"))) stop ("Invalid 'pred.type' argument.")
if (pred.type == "F" && family(model)$family != "binomial") stop ("pred.type='F' is only applicable to binomial models.")
vars <- names(model$coefficients)[-1]
if (is.null(groups)) stop ("If you provide a 'model' argument, you must also provide 'groups'.")
if (!all(vars %in% groups[ , 1])) stop ("All variables in 'model' must exist (with the same exact name) in the first column of 'groups'.")
if (length(unique(groups[ , 2])) != length(unique(trimws(groups[ , 2])))) warning ("Some values in 'groups' have leading or trailing spaces and are treated separately; consider using trimws() first.")
groups <- groups[groups[ , 1] %in% vars, ]
if (length(unique(groups[ , 2])) > 3) stop ("All variables in 'model' must belong to a maximum of three different 'groups'.")
factors <- unique(groups[ , 2])
factors_ABC <- LETTERS[1:length(factors)]
groups[ , 3] <- factors_ABC[match(groups[ , 2], factors)]
c1 <- factors_ABC # each individual factor
if (length(factors) > 1) c2 <- as.data.frame(combn(factors_ABC, 2)) else c2 <- NA # each combination of 2 factors
if (length(factors) > 2) c3 <- as.data.frame(combn(factors_ABC, 3)) else c3 <- NA # each combination of 3 factors
combinations <- unname(c(c1, c2, c3))
combinations <- Filter(Negate(anyNA), combinations) # like na.omit for list
mods <- vector("list", length(combinations))
names(mods) <- sapply(combinations, paste, collapse = "")
response <- names(model$model)[1]
for (fac in 1:length(combinations)) {
vars_fac <- groups[groups[ , 3] %in% unlist(combinations[fac]), 1]
form <- as.formula(paste(response, "~", paste(vars_fac, collapse = "+")))
mods[[fac]] <- glm(form, data = model$model, family = family(model)$family)
}
preds <- data.frame(matrix(nrow = nrow(model$model), ncol = 0))
if (pred.type == "Y") type <- "link" else type <- "response"
for (m in names(mods)) {
preds[ , m] <- predict(mods[[m]], type = type)
}
if (pred.type == "F") {
n1 <- sum(model$y == 1)
n0 <- sum(model$y == 0)
preds <- as.data.frame(sapply(preds, function(p) (p/(1-p)) / ((n1/n0) + (p/(1-p))))) # can't import from fuzzySim because fuzzySim imports modEvA
}
rsq <- rep(NA_real_, length(mods))
names(rsq) <- names(mods)
for (f in names(preds)) {
# if (family(model)$family == "gaussian") linmod <- lm(model$model[ , 1] ~ preds[ , f])
# else linmod <- lm(preds[ , ncol(preds)] ~ preds[ , f])
# rsq[f] <- summary(linmod)$r.squared
if (family(model)$family == "gaussian") corr <- cor(model$model[ , 1], preds[ , f], method = cor.method)
else corr <- cor(preds[ , ncol(preds)], preds[ , f], method = cor.method)
rsq[f] <- corr ^ 2
}
}
if (is.null(model)) {
partials <- c(A, B, C, AB, BC, AC, ABC)
if (!all(na.omit(partials) >= 0 & na.omit(partials) <= 1)) stop ("Values of A, B, etc. must be between 0 and 1.")
if (all(is.finite(partials[c(1:2, 4)])) && all(is.na(partials[c(3, 5:7)])))
twofactors <- TRUE
else if (all(is.finite(partials)))
twofactors <- FALSE
else stop ("You must provide numeric values for either A, B and AB (for variation partitioning among two factors) or A, B, C, AB, BC, AC and ABC (for variation partitioning among three factors). See Details.")
} else {
if (length(rsq) <= 4) twofactors <- TRUE else twofactors <- FALSE
A <- rsq["A"]
B <- rsq["B"]
C <- rsq["C"]
AB <- rsq["AB"]
AC <- rsq["AC"]
BC <- rsq["BC"]
ABC <- rsq["ABC"]
A.name <- unique(groups[groups[ , 3] == "A", 2])
B.name <- unique(groups[groups[ , 3] == "B", 2])
C.name <- unique(groups[groups[ , 3] == "C", 2])
}
totalexpl <- ifelse(twofactors, AB, ABC)
unexpl <- 1 - totalexpl
if (twofactors) {
Apure <- totalexpl - B
Bpure <- totalexpl - A
ABoverlap <- totalexpl - Apure - Bpure
output.names <- c(A.name, B.name, paste(A.name, B.name, sep = "_"),
"Unexplained")
results <- data.frame(c(Apure, Bpure, ABoverlap, unexpl),
row.names = output.names)
} else { # end if 2 factors
Apure <- totalexpl - BC
Bpure <- totalexpl - AC
Cpure <- totalexpl - AB
ABoverlap <- totalexpl - Apure - Bpure - C
BCoverlap <- totalexpl - Bpure - Cpure - A
ACoverlap <- totalexpl - Apure - Cpure - B
ABCoverlap <- totalexpl - Apure - Bpure - Cpure - ABoverlap - BCoverlap - ACoverlap
output.names <- c(paste(A.name),
paste(B.name),
paste(C.name),
paste(A.name, B.name, sep = "_"),
paste(B.name, C.name, sep = "_"),
paste(A.name, C.name, sep = "_"),
paste(A.name, B.name, C.name, sep = "_"),
"Unexplained")
results <- data.frame(c(Apure, Bpure, Cpure, ABoverlap, BCoverlap,
ACoverlap, ABCoverlap, unexpl),
row.names = output.names)
} # end else
colnames(results) <- "Proportion"
#n <- nrow(results)
#if(model.type == "GLM") results <- results[1:(n-1),] # deletes "unexplained" line (data unavailable for GLM)
if (plot) { # adapted from Daniel's http://stackoverflow.com/questions/1428946/venn-diagrams-with-r
circle <- function(x, y, r, col = NA) {
ang <- seq(0, 2 * pi, length = 100)
xx <- x + r * cos(ang)
yy <- y + r * sin(ang)
polygon(xx, yy, col = col)
} # end circle funtion (by Daniel)
Apure <- round(Apure, plot.digits) # shorten values for plotting
Bpure <- round(Bpure, plot.digits)
ABoverlap <- round(ABoverlap, plot.digits)
if(!twofactors) {
Cpure <- round(Cpure, plot.digits)
BCoverlap <- round(BCoverlap, plot.digits)
ACoverlap <- round(ACoverlap, plot.digits)
ABCoverlap <- round(ABCoverlap, plot.digits)
}
if (twofactors) {
plot(0, 0, ylim = c(-1, 10), xlim = c(-1, 10), type = "n", axes = FALSE,
ylab = "", xlab = "", main = main, cex.main = cex.main)
circle(4.5, 3, 3, col = ifelse(colr, rgb(1, 0, 0, 0.5), NA))
circle(4.5, 6, 3, col = ifelse(colr, rgb(0, 1, 0, 0.5), NA))
text(x = c(4.5, 4.5), y = c(9.5, -0.5), labels = c(A.name, B.name),
cex = cex.names)
text(x = c(4.5, 4.5, 4.5), y = c(7, 4.75, 2), c(Apure, ABoverlap, Bpure),
cex = cex.values)
} else { # end if 2 factors
plot(0, 0, ylim = c(-1, 10), xlim = c(-1, 10), type = "n", axes = FALSE,
ylab = "", xlab = "", main = main, cex.main = cex.main)
circle(3, 6, 3, col = ifelse(colr, rgb(1, 0, 0, 0.5), NA))
circle(6, 6, 3, col = ifelse(colr, rgb(0, 1, 0, 0.5), NA))
circle(4.5, 3, 3, col = ifelse(colr, rgb(0, 0, 1, 0.5), NA))
#Cname.loc = ifelse((plot.unexpl), 6, 4.5)
text(x = c(2.5, 6.5, 4.5), y = c(9.5, 9.5, -0.5),
labels = c(A.name, B.name, C.name), cex = cex.names, adj = c(0.5, 0.5, 0))
text(x = c(1.8, 7.2, 4.5, 4.5, 2.8, 6.2, 4.5), y = c(6.6, 6.6, 2, 7, 4, 4, 5), labels = c(Apure, Bpure, Cpure, ABoverlap, ACoverlap, BCoverlap, ABCoverlap), cex = cex.values)
} # end if 2 factors else
if (plot.unexpl) {
rect(-1, -1, 10, 10)
text(x = -0.9, y = -0.2, label = paste0("Unexplained\n", round(unexpl, plot.digits)), adj = 0, cex = cex.values)
}
} # end if plot
if (all.equal(sum(results, na.rm = TRUE), 1)) cat("")
else warning ("Results don't sum up to 1; are you sure your input data are correct?") # but this doesn't work because results always sum to 1 anyway
if (!return.models) return(results)
return(list(models = mods, varpart = 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.