# Rbrul version 2.3 # Daniel Ezra Johnson # Oct 4 2015 w <- options("warn") options(warn = 2) if (.Platform$OS.type == "unix") a <- try(system("open http://www.danielezrajohnson.com/analytics.html", ignore.stderr = T, wait = F), silent = T) else a <- try(system("start http://www.danielezrajohnson.com/analytics.html", ignore.stderr = T, wait = F, invisible = T, show.output.on.console = F), silent = T) rm("a", envir = globalenv()) options(warn = as.numeric(w)) rbrul <- function(packages = T){ w <- options("warn") options(warn = 2) pdf <- print.data.frame unlockBinding("print.data.frame", asNamespace("base")) assign("print.data.frame", rownames.off, envir = asNamespace("base")) cntrsts <- options("contrasts") rb <- try(rbrul2(packages = packages)) if ("try-error" %in% class(rb)){ cat(sep = "\n", "\nSorry, there was an error and Rbrul crashed.") cat(sep = "\n", "Help is available by emailing danielezrajohnson@gmail.com.")} assign("print.data.frame", pdf, envir = asNamespace("base")) lockBinding("print.data.frame", asNamespace("base")) options(contrasts = cntrsts[[1]]) options(warn = as.numeric(w)) } rbrul2 <- function(packages){ cat("\n", "Welcome to Rbrul version 2.3 - Cleveland Circle") si <- try(sessionInfo(), silent = T) cat("\n", "Your current R version is", paste(si$R.version$major, ".", si$R.version$minor, " -", sep=""), si$R.version$nickname) if (as.numeric(si$R.version$major) < 3){ cat(sep = "\n", "Please go to http://cran.r-project.org/ and update to the latest version of R.") return()} cat("\n") cat("\n", "In case of trouble, make sure your R version is updated, and try loading Rbrul again.") cat("\n", "Please don't hesitate to email danielezrajohnson@gmail.com for help. (This works.)") cat("\n") cat("\n", "Windows users may need to run R as administrator the first time they use Rbrul,") cat("\n", "or when updates occur. Right-click to do this. You may need an admin password.") cat("\n") cat("\n", "If you get an error saying 'invalid multibyte string', it means foreign characters") cat("\n", "are in your data: accents, diacritics, digraphs, or even worse! Get rid of them!") cat("\n") cat("\n", "Rbrul is no longer suppressing (red) convergence warnings. Models that didn't") cat("\n", "converge properly may be too complex for the data. Trust them at your own risk.") cat("\n") get.packages(packages) vars <<- c("addmargins", "adjusting.pick", "apps", "blups", "centered", "combo", "continuous", "continuous.pick", "contr", "datafile", "decimals", "down", "down.model", "down.preds", "fam", "exit.pick", "familia", "fixed", "fullpath", "gold", "interactions", "interactions.picks", "loaded.flag", "main.pick", "model", "modeling.pick", "nloptr.default.options", "non.apps", "nopackages", "old.path", "old.predictors", "old.response", "one", "one.model", "one.preds", "p.flag", "predictors", "predictors.pick", "proportion", "random.i", "random.s", "response", "response.pick", "sepchar", "silent", "thresh", "threshold", "up", "up.model", "up.preds", "variables", "verbose", "w") vars <<- c(vars, c("adjusting", "bin2num", "crosstab", "count", "cross", "crosstabs", "custom.plot", "data.structure", "datafile2", "exclude", "extract", "get.contrasts", "get.continuous", "get.data", "get.input", "get.interactions", "get.p", "get.packages", "get.predictors", "get.random", "get.response", "get.settings", "get.variables", "fit.model", "get.menu", "get.weights", "load.variables", "maxfun", "modeling", "panel.rsuperpose", "panel.rxyplot", "partial", "plotting", "plotting.pick", "rbrul", "rbrul2", "reclass", "recode", "relev", "rename", "reset", "retain", "rownames.off", "save.data", "set.variables", "shortpath", "split.term", "star", "step.down", "step.one", "step.up", "step.updown", "sumify", "test.chisq", "tolerance", "transform", "trim", "vif.dej")) for (v in vars) overwrite(v, NULL) overwrite("silent", FALSE) overwrite("verbose", T) overwrite("gold", T) overwrite("centered", T) overwrite("blups", T) overwrite("thresh", 0.05) overwrite("decimals", 3) overwrite("maxfun", 100000) overwrite("contr", "contr.sum") options(contrasts = c(unordered = contr, ordered = "contr.poly")) overwrite("loaded.flag", F) overwrite("sepchar", ", ") overwrite("fullpath", "") overwrite("shortpath", paste(getwd(), "/", sep = "")) overwrite("tolerance", 0.002) old.path <<- shortpath main.pick <<- 99 while (main.pick != 0){ cat("\n") if (is.null(datafile)) { cat("\n") cat(sep = "\n", "No data loaded.") cat("\n")} else { cat("\n") cat(sep = "\n", paste("Current data file is:", fullpath)) data.structure()} cat(sep = "\n", "MAIN MENU") if (!is.null(datafile)) cat(sep = "\n", "1-load/save data 2-adjust data") else cat(sep = "\n", "1-load/save data") if (!is.null(datafile)) cat(sep = "\n", "4-crosstabs 5-modeling 6-plotting") if (!is.null(datafile)) cat(sep = "\n", "8-restore data 9-reset 0-exit") else cat(sep = "\n", "9-reset 0-exit") main.pick <<- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(main.pick) == 0) main.pick <<- 99 if (main.pick == 1) get.data() if (main.pick == 2) if (!is.null(datafile)) adjusting() else cat(sep = "\n", "Please load data first.") if (main.pick == 4) if (!is.null(datafile)) crosstabs() else cat(sep = "\n", "Please load data first.") if (main.pick == 5) if (!is.null(datafile)) modeling() else cat(sep = "\n", "Please load data first.") if (main.pick == 6) if (!is.null(datafile)) plotting() else cat(sep = "\n", "Please load data first.") if (main.pick == 8) if (!is.null(datafile)) {datafile <<- datafile2; loaded.flag <<- T} else cat(sep = "\n", "No data loaded!") if (main.pick == 9) reset()} cat(sep = "\n", "Discard Rbrul from the workspace? (d-discard, Enter-keep)") repeat {exit.pick <- try(as.character(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(exit.pick) == 0){ setwd(old.path) break} if (exit.pick %in% c("d", "D")) { for (v in vars) try(rm(list = eval(v), envir = globalenv()), silent = T) try(rm("vars", envir = globalenv()), silent = T) try(rm("overwrite", envir = globalenv()), silent = T) # re sorce if reloaded? break} cat(sep = "\n", "Invalid entry, please try again.")}} addmargins <- function (A, margin = 1:length(dim(A)), FUN = sum, quiet = FALSE) { if (is.null(dim(A))) stop("'A' must be an array or table") n.sid <- length(margin) miss.FUN <- missing(FUN) if (length(FUN) == 1 && !is.list(FUN)) { fname <- if (!miss.FUN) deparse(substitute(FUN)) else "total" FUN <- list(FUN) names(FUN) <- fname } if (!miss.FUN) { add.names <- function(thelist) { n <- names(thelist) if (is.null(n)) n <- rep("", length(thelist)) for (i in seq_along(thelist)[-1]) { if (!is.call(thelist[[i]])) { if (n[i] == "") n[i] <- as.character(thelist[[i]]) } else if (as.character(thelist[[i]][[1]]) == "list") thelist[[i]] <- add.names(thelist[[i]]) } names(thelist) <- n thelist } if (mode(substitute(FUN)) == "call") FUN <- eval(add.names(substitute(FUN))) if (is.null(names(FUN))) names(FUN) <- rep("", length(FUN)) } if (length(FUN) != n.sid) { if (length(FUN) == 1) FUN <- rep(FUN, n.sid) else stop(gettextf("length of FUN, %d, \n does not match the length of the margins, %d", length(FUN), n.sid), domain = NA) } fnames <- vector("list", n.sid) for (i in seq_along(FUN)) { fnames[[i]] <- names(FUN)[i] if (is.list(FUN[[i]])) { topname <- fnames[[i]] fnames[[i]] <- names(FUN[[i]]) blank <- fnames[[i]] == "" fnames[[i]][blank] <- seq_along(blank)[blank] if (topname == "") { fnames[[i]][blank] <- paste("Margin ", margin[i], ".", fnames[[i]][blank], sep = "") } else { fnames[[i]] <- paste(topname, ".", fnames[[i]], sep = "") } } else if (fnames[[i]] == "") fnames[[i]] <- paste("Margin", margin[i]) } expand.one <- function(A, margin, FUN, fnames) { if (class(FUN) != "list") FUN <- list(FUN) d <- dim(A) n.dim <- length(d) n.mar <- length(FUN) newdim <- d newdim[margin] <- newdim[margin] + n.mar newdimnames <- dimnames(A) newdimnames[[margin]] <- c(newdimnames[[margin]], fnames) n.new <- prod(newdim) skip <- prod(d[1:margin]) runl <- skip/d[margin] apos <- rep(c(rep(TRUE, skip), rep(FALSE, n.mar * runl)), n.new/(skip + n.mar * runl)) values <- numeric(length(apos)) values[apos] <- as.vector(A) for (i in 1:n.mar) { mtab <- if (n.dim > 1) { apply(A, (1:n.dim)[-margin], FUN[[i]]) } else FUN[[i]](A) select <- rep(FALSE, n.mar) select[i] <- TRUE mpos <- rep(c(rep(FALSE, skip), rep(select, each = runl)), prod(dim(A))/skip) values[mpos] <- as.vector(mtab) } new.A <- array(values, dim = newdim, dimnames = newdimnames) if (inherits(A, "table")) class(new.A) <- c("table", class(new.A)) new.A } new.A <- A for (i in 1:n.sid) new.A <- expand.one(A = new.A, margin = margin[i], FUN = FUN[[i]], fnames = fnames[[i]]) if (!quiet && !miss.FUN && n.sid > 1) { cat("Margins computed over dimensions\nin the following order:\n") for (i in 1:n.sid) cat(paste(i), ": ", names(dimnames(A))[margin[i]], "\n", sep = "") } new.A } adjusting <- function(){ adjusting.pick <<- 99 while (!adjusting.pick %in% c(9, 0)){ cat("\n") cat(sep = "\n", "ADJUSTING MENU") cat(sep = "\n", "1-change class 2-rename 3-exclude 4-retain 5-recode") cat(sep = "\n", "6-relevel 7-center/transform 8-count 9-main menu 0-exit") cat(sep = "\n", "10-make full interaction group 11-make partial interaction group") adjusting.pick <<- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(adjusting.pick) == 0) adjusting() if (adjusting.pick == 1) reclass() if (adjusting.pick == 2) rename() if (adjusting.pick == 3) exclude() if (adjusting.pick == 4) retain() if (adjusting.pick == 5) recode() if (adjusting.pick == 6) relev() if (adjusting.pick == 7) transform() if (adjusting.pick == 8) count() if (adjusting.pick == 10) star() if (adjusting.pick == 11) partial() if (adjusting.pick == 0) main.pick <<- 0}} bin2num <- function(x){ if (class(x)!= "factor") return(x) if (length(levels(x))!= 2) return(x) levels(x)[1] <- 0 levels(x)[2] <- 1 as.numeric(as.character(x))} count <- function(){ cat(sep = "\n", paste("Column(s) to count? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {count.pick <- try(as.numeric(scan(what = "character", nmax = length(names(datafile)), quiet = T)), silent = T) countees <- names(datafile)[count.pick] if (!NA %in% countees) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(countees) > 0) for (i in countees){ count.col <- paste("count", i, sep = ".") for (j in unique(as.character(datafile[, i]))) datafile[datafile[, i] == j, count.col] <<- length(datafile[datafile[, i] == j, 1])}} cross <- function(quiet = F, crossees = c(NULL, NULL)){ if (!quiet){ cat(sep = "\n", paste("Choose two variables to cross? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {cross.pick <- try(as.numeric(scan(what = "character", nmax = 2, quiet = T)), silent = T) crossees <- names(datafile)[cross.pick] if (length(crossees) == 2 & !NA %in% crossees) break cat(sep = "\n", "Invalid entry, please try again.")}} x <- datafile[, crossees[1]] y <- datafile[, crossees[2]] z <- list("error") zz <- vector() if (class(x) == "character") x <- as.factor(x) if (class(y) == "character") y <- as.factor(y) if (class(x) == "factor" & class(y) == "factor"){ x <- x[drop = T] y <- y[drop = T] newfac <- 1 for (i in 1:length(contrasts(x)[1, ])){ for (j in 1:length(contrasts(y)[1, ])){ z[[newfac]] <- get.contrasts(x, i)*get.contrasts(y, j) newfac <- newfac + 1}}} if (class(x) == "factor" & class(y) %in% c("integer", "numeric")| class(y) == "factor" & class(x) %in% c("integer", "numeric")){ if (class(y) == "factor") {t <- y; y <- x; x <- t} for (i in 1:length(contrasts(x)[1, ])) z[[i]] <- get.contrasts(x, i)*y} if (class(x) %in% c("integer", "numeric") & class(y) %in% c("integer", "numeric")) z <- list(x*y) for (i in 1:length(z)){ if (length(z) == 1) datafile[, zz[i] <-paste(crossees[1], ".", crossees[2], sep = "")] <<- z[[i]] else datafile[, zz[i] <-paste(crossees[1], ".", crossees[2], ".", letters[i], sep = "")] <<- z[[i]]} zz} crosstabs <- function(){ cat(sep = "\n", paste("Cross-tab: factor for columns?", get.menu(names(datafile)))) repeat {columns.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) columns <- names(datafile)[columns.pick] if (!NA %in% columns & length(columns.pick) > 0) { datafile[, columns] <- as.factor(datafile[, columns]); break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", paste("Cross-tab: factor for rows? (", get.menu(names(datafile)[!names(datafile) %in% columns], which(!names(datafile) %in% columns), paren = F), " Enter-none)", sep = "")) repeat {rows.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) rows <- names(datafile)[rows.pick] if (length(rows.pick) == 0) {rows <- 1; break} if (!NA %in% rows) { datafile[, rows] <- as.factor(datafile[, rows]); break} cat(sep = "\n", "Invalid entry, please try again.")} if (rows == 1) pages <- 1 else { cat(sep = "\n", paste("Cross-tab: factor for 'pages'? (", get.menu(names(datafile)[!names(datafile) %in% c(columns, rows)], which(!names(datafile) %in% c(columns, rows)), paren = F), " Enter-none)", sep = "")) repeat {pages.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) pages <- names(datafile)[pages.pick] if (length(pages.pick) == 0) {pages <- 1; break} if (!NA %in% pages) { datafile[, pages] <- as.factor(datafile[, pages]); break} cat(sep = "\n", "Invalid entry, please try again.")}} if (!is.null(response)){ if (response %in% names(datafile)){ cat(sep = "\n", paste("Cross-tab: variable for cells? (1-proportion/mean of ", response, " 2-change response variable Enter-counts)", sep = ""))} else{ cat(sep = "\n", "Cross-tab: variable for cells? (1-response proportion/mean Enter-counts)")}} else{ cat(sep = "\n", "Cross-tab: variable for cells? (1-response proportion/mean Enter-counts)")} repeat {cells.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(cells.pick) == 0){ cat(sep = "\n", "counts") crosstab <<- addmargins(xtabs(as.formula(paste("~", rows, "+", columns, "+", pages, sep = "")), datafile)) return(print(crosstab))} if (cells.pick == 1){ if (is.null(response)) get.response() else if (!response %in% names(datafile)) get.response() break} if (cells.pick == 2){ get.response() break} cat(sep = "\n", "Invalid entry, please try again.")} if (class(datafile[, response]) %in% c("integer", "numeric")) {cat(sep = "\n", paste("mean of", response)) crosstab <<- round(addmargins(xtabs(as.formula(paste(response, "~", rows, "+", columns, "+", pages, sep = "")), datafile))/ addmargins(xtabs(as.formula(paste("~", rows, "+", columns, "+", pages, sep = "")), datafile)), decimals) return(print(crosstab))} else {if (is.null(apps)) get.response() cat(sep = "\n", paste("proportion of", response, " = ", apps)) crosstab <<- round(addmargins(xtabs(as.formula(paste("~", rows, "+", columns, "+", pages, sep = "")), datafile[datafile[, response] %in% apps, ]))/ addmargins(xtabs(as.formula(paste("~", rows, "+", columns, "+", pages)), datafile)), decimals)} return(print(crosstab))} custom.plot <- function(){ pal <- F data.structure() cat(sep = "\n", paste("Data variables:", get.menu(names(datafile), paren = F))) cat("\n") if (is.null(model)){ cat(sep = "\n", "No model loaded."); combo <- datafile} else { mod <- model[[1]] if ("glm" %in% class(mod)){ modelvars <- mod$data modelvars$predictions <- mod$fitted.values if (mod$family[1] == "gaussian"){ familia <- "linear" modelvars$residuals <- mod$residuals} else familia <- "logistic"} if (class(mod)[1] %in% c("lmerMod", "glmerMod")){ modelvars <- mod@frame if (isGLMM(mod) == F){ modelvars$residuals <- residuals(mod) modelvars$predictions <- getME(mod, "X") %*% fixef(mod) familia <- "linear"} else {modelvars$predictions <- plogis(getME(mod, "X") %*% fixef(mod)) familia <- "logistic"} modelvars$predictions.including.random.effects <- fitted(mod)} cat(sep = "\n", paste("Model variables:", get.menu(names(modelvars), index = (length(names(datafile))+1):(length(names(datafile))+length(names(modelvars))), paren = F))) if (length(datafile[, 1]) == length(modelvars[, 1])) combo <- cbind(datafile, modelvars) else { cat("\n") cat(sep = "\n", "Warning! Data file and model do not match. Use model variables only.") combo <- modelvars for (i in 1:length(names(datafile))) combo <- cbind(dummy = 1, combo)}} n <- length(combo[, 1]) cat("\n") cat(sep = "\n", "Choose variable for y-axis?") repeat{y.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) y <- names(combo)[y.pick] if (length(y) == 0) next if (!NA %in% y & y != "dummy") break cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Choose variable for x-axis?") repeat{x.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) x <- names(combo)[x.pick] if (length(x) == 0) next if (!NA %in% x & x != "dummy") break cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Separate (and color) data according to the values of which variable? (press Enter to skip)") repeat{g.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) g <- names(combo)[g.pick] if (length(g) == 0) {g <- NULL; gmean.pick <- 0; break} if (!NA %in% g & g != "dummy") { cat(sep = "\n", paste("Also show data (in black) averaged over all values of ", g, "? (1-yes Enter-no)", sep = "")) repeat {gmean.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(gmean.pick) == 0){ gmean.pick <- 0 break} if (gmean.pick == 1){ combo2 <- combo combo2[, g] <- "XXX" combo <- rbind(combo, combo2) pal <- T break} cat(sep = "\n", "Invalid entry, please try again.")} break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Split data into horizontal panels according to which variable? (press Enter to skip)") repeat{cols.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) cols <- names(combo)[cols.pick] if (length(cols) == 0) {cols <- 1; break} if (!NA %in% cols & cols!= "dummy") break cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Split data into vertical panels according to which variable? (press Enter to skip)") repeat{rows.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) rows <- names(combo)[rows.pick] if (length(rows) == 0) {rows <- 1; break} if (!NA %in% rows & rows!= "dummy") break cat(sep = "\n", "Invalid entry, please try again.")} if (cols == 1 & rows == 1) strip <- F else strip <- T cat(sep = "\n", "Type of points to plot? (raw points not recommended for binary data)", "(0-no points 1-raw points Enter-mean points)") repeat{points.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(points.pick) == 0) {points.pick <- 2; p <- "pm"; break} if (points.pick == 0) {p <- NULL; break} if (points.pick == 1) {p <- "p"; break} cat(sep = "\n", "Invalid entry, please try again.")} scale.factor <- 0 if ((gmean.pick == 1 & points.pick == 1) | points.pick == 2){ cat(sep = "\n", "Scale points according to the number of observations?", "Enter size factor between 0.1 and 10 (1 = Enter = default)", "or 0 to not scale points.") repeat{scale.factor <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(scale.factor) == 0) {scale.factor <- 1; break} if (scale.factor == 0) break if (0.1 <= scale.factor & scale.factor <= 10) break cat(sep = "\n", "Invalid entry, please try again.")}} cat(sep = "\n", "Type of lines to plot (raw lines not recommended for binary data)?", "0-no lines 1-raw lines Enter-mean lines)") repeat{lines.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(lines.pick) == 0) {l <- "lm"; break} if (lines.pick == 0) {l <- NULL; break} if (lines.pick == 1) {l <- "l"; break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Add a reference line? (1-diagonal [y = x] 2-horizontal [y = 0] Enter-none)") repeat{refline.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(refline.pick) == 0) {refline <- ""; break} if (refline.pick == 1) {refline <- "d"; break} if (refline.pick == 2) {refline <- "h"; break} cat(sep = "\n", "Invalid entry, please try again.")} combo[, y] <- bin2num(combo[, y]) if (length(g) == 0) gg <- "NULL" else gg <- paste("combo$", g, sep = "") form <- as.formula(paste(y, "~", x, "|", cols, "+", rows, sep = "")) if (strip) if (pal){ superpose.symbol.old <- trellis.par.get("superpose.symbol") superpose.line.old <- trellis.par.get("superpose.line") superpose.symbol <- superpose.symbol.old superpose.line <- superpose.line.old superpose.line$col <- c(rainbow(length(unique(combo[, g]))-1), "black") superpose.symbol$col <- c(rainbow(length(unique(combo[, g]))-1), "black") trellis.par.set(superpose.symbol = superpose.symbol) trellis.par.set(superpose.line = superpose.line)} go.flag <- T if (!is.null(familia)) { if (familia == "logistic" & length(grep("predictions", x)) > 0){ go.flag <- F print(xyplot(form, groups = eval(parse(text = gg)), type = c(l, p), auto.key = T, data = combo, xlim = c(-0.07, 1.07), panel = function(...){ if (refline == "d") panel.abline(0, 1) if (refline == "h") panel.abline(h = 0) panel.rxyplot(...)}, strip = strip, scale.factor = scale.factor, scales = list(alternating = 1)))}} if (go.flag){ print(xyplot(form, groups = eval(parse(text = gg)), type = c(l, p), auto.key = T, data = combo, panel = function(...){ if (refline == "d") panel.abline(0, 1) if (refline == "h") panel.abline(h = 0) panel.rxyplot(...)}, strip = strip, scale.factor = scale.factor, scales = list(alternating = 1)))} if (pal){ trellis.par.set(superpose.symbol = superpose.symbol.old) trellis.par.set(superpose.line = superpose.line.old)}} # layout for rows and columns # automatically order factors like in B&P # prompt for y = x line if y and x are in predicted, response data.structure <- function(){ cat("\n") cat(sep = "\n", "Current data structure:") classes <- sapply(datafile, class) for (i in 1:length(classes)){ cat(paste(names(datafile)[i], " (", classes[i], " with ", l <- length(u <- unique(datafile[, i])), ifelse(l == 1, " value): ", " values): "), sep = "")) cat(sep = "\n", paste(paste(u[1:ifelse(l < 5, l, 5)], collapse = " "), ifelse(l < 6, "", "..."), collapse = ""))} cat(sep = "\n", paste("Total tokens:", length(datafile[, 1]), sep = " ")) cat("\n")} exclude <- function(){ cat(sep = "\n", paste("Predictor to exclude from?", get.menu(names(datafile)), sep = "")) repeat {exclude.group.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) exclude.group <- names(datafile)[exclude.group.pick] if (!NA %in% exclude.group) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(exclude.group) == 0) return() if (class(datafile[, exclude.group]) %in% c("factor", "character")){ fac <- levels(as.factor(datafile[, exclude.group])) cat(sep = "\n", paste("Exclude which factors from ", exclude.group, "?", get.menu(fac), sep = "")) repeat {exclude.factors.pick <- try(as.numeric(scan(what = "character", nmax = length(fac), quiet = T)), silent = T) exclude.factors <- fac[exclude.factors.pick] if (!NA %in% exclude.factors) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(exclude.factors) == 0) return() datafile <<- subset(datafile, !datafile[, exclude.group] %in% exclude.factors) datafile[, exclude.group] <<- datafile[, exclude.group][drop = T]} if (class(datafile[, exclude.group]) %in% c("numeric", "integer")){ r <- range(datafile[, exclude.group]) cat(sep = "\n", paste("Range of ", exclude.group, " is ", r[1], " to ", r[2], ". Choose cut-off value.", sep = "")) repeat {exclude.cutoff <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (!NA %in% exclude.cutoff) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(exclude.cutoff) == 0) return() cat(sep = "\n", paste("Exclude values of ", exclude.group, " that are a-less than ", exclude.cutoff, ", b-greater than ", exclude.cutoff, ", c-equal to ", exclude.cutoff, "?", sep = "")) repeat {exclude.dir <- try(as.character(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (exclude.dir %in% c("a", "b", "c", "A", "B", "C")) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(exclude.dir) == 0) return() if (exclude.dir %in% c("a", "A")) datafile <<- subset(datafile, datafile[, exclude.group] >= exclude.cutoff) if (exclude.dir %in% c("b", "B")) datafile <<- subset(datafile, datafile[, exclude.group] <= exclude.cutoff) if (exclude.dir %in% c("c", "C")) datafile <<- subset(datafile, datafile[, exclude.group] != exclude.cutoff)} datafile <<- droplevels(datafile)} extract <- function(model){ if ("try-error" %in% class(model)) return(cat("Error Message:\n", model, "\n")) if ("list" %in% class(model)) { cat(sep = "\n", model[[2]]); cat("\n") model <- model[[1]]} extraction <- list() group.labels <- vector() factor.labels <- vector() contrastss <- vector() groups <- vector() if (is.null(contr)) contr <<- "contr.treatment" n <- length(fitted(model)) if ("glm" %in% class(model)){ term.labels <- attr(model$terms, "term.labels") factor.slot <- "model$model$" intercept <- unname(model$coefficients[1]) coefficientss <- model$coefficients[-1]} if (any(c("lmerMod", "glmerMod") %in% class(model)[1])){ term.labels <- attr(attr(model@frame, "terms"), "term.labels") term.labels <- term.labels[!term.labels %in% names(VarCorr(model))] factor.slot <- "model@frame$" contrr <- unname(unlist(attr(getME(model, "X"), "contrasts")[1])) intercept <- unname(fixef(model)[1]) coefficientss <- fixef(model)[-1]} if (length(term.labels) > 0){ for (t in term.labels){ if (all(split.term(t) %in% continuous)){ factor.labels <- c(factor.labels, "+1") group.labels <- c(group.labels, t) contrastss <- c(contrastss, coefficientss[1]) coefficientss <- coefficientss[-1]} if (all(!split.term(t) %in% continuous)){ num.coefs <- vector() name.coefs <- vector() names.coefs <- data.frame() for (f in split.term(t)){ name.coefs <- levels(as.factor(eval(parse(text = paste(factor.slot, f, sep = ""))))) num.coefs <- c(num.coefs, length(name.coefs)-1) names.coefs <- rbind(names.coefs, cbind(factorr = f, level = name.coefs, first = c(T, rep(F, length(name.coefs)-1)), last = c(rep(F, length(name.coefs)-1), T)))} factor.update <- levels(as.factor(eval(parse(text = paste(factor.slot, split.term(t), sep = "", collapse = ":"))))) factor.labels <- c(factor.labels, factor.update) group.labels <- c(group.labels, rep(t, length(factor.update))) if (contr == "contr.treatment"){ for (f in factor.update){ base.levels <- names.coefs$level[names.coefs$first == T] if (any(split.term(f) %in% base.levels)) contrastss <- c(contrastss, 0) else { contrastss <- c(contrastss, coefficientss[1]) coefficientss <- coefficientss[-1]}}} if (contr == "contr.sum"){ contrastsss <- vector() for (f in factor.update){ omitted.levels <- names.coefs$level[names.coefs$last == T] if (any(split.term(f) %in% omitted.levels)) contrastsss <- c(contrastsss, NA) else { contrastsss <- c(contrastsss, coefficientss[1]) coefficientss <- coefficientss[-1]}} contrastss <- c(contrastss, sumify(contrastsss, length(levels(as.factor(eval(parse(text = paste(factor.slot, split.term(t)[1])))))))) }} if (!(all(split.term(t) %in% continuous)) & !(all(!split.term(t) %in% continuous))){ factor.update <- levels(as.factor(eval(parse(text = paste( factor.slot, split.term(t)[!split.term(t) %in% continuous], sep = ""))))) if (split.term(t)[1] %in% continuous) factor.labels <- c(factor.labels, paste("+1", factor.update, sep = ":")) else factor.labels <- c(factor.labels, paste(factor.update, "+1", sep = ":")) group.labels <- c(group.labels, rep(t, length(factor.update))) if (contr == "contr.treatment") contrastss <- c(contrastss, 0, coefficientss[1:length(factor.update)-1]) if (contr == "contr.sum") contrastss <- c(contrastss, s <-coefficientss[1:length(factor.update)-1], -sum(s)) coefficientss <- coefficientss[-c(1:length(factor.update)-1)]}} fram <- data.frame(group = group.labels, factor = factor.labels, coefs = contrastss) groups <- unique(as.character(fram$group)) for (i in 1:(length(groups))){ extraction[[i]] <- fram[fram$group == groups[i], -1] names(extraction)[i] <- groups[i] names(extraction[[i]]) <- c("factor", "coefs") for (j in 1:length(extraction[[i]]$factor)){ if ("/" %in% split.term(as.character(extraction[[i]]$factor[j]))){ extraction[[i]]$coefs[j] <- NA}} if (groups[i] %in% continuous) names(extraction[[i]]) <- c("continuous", "coefs") if (length(split.term(groups[i])) > 1){ names(extraction[[i]]) <- c("factor:factor", "coefs") if (all(split.term(groups[i]) %in% continuous)) names(extraction[[i]]) <- c("continuous:continuous", "coefs") else { if (split.term(groups[i])[1] %in% continuous) names(extraction[[i]]) <- c("continuous:factor", "coefs") if (split.term(groups[i])[2] %in% continuous) names(extraction[[i]]) <- c("factor:continuous", "coefs")}} if (names(extraction[[i]])[1] %in% c("factor", "factor:factor")){ for (j in 1:length(extraction[[i]][, 1])){ extraction[[i]][j, "tokens"] <- length(which(as.character(eval(parse(text=paste(factor.slot, split.term(groups[i]),sep="",collapse=":")))) == as.character(extraction[[i]][j, 1]))) if (fam == "binomial") extraction[[i]][j,proportion] <- round(length(which(as.character(eval(parse(text=paste( factor.slot,split.term(groups[i]),sep="",collapse=":")))) == as.character(extraction[[i]][j, 1]) & as.character(eval(parse(text = paste(factor.slot, response, sep = "")))) == apps)) / extraction[[i]][j, "tokens"], decimals) if (fam == "gaussian") if (length(split.term(groups[i])) == 1) extraction[[i]][j, proportion] <- round(mean(datafile[, response][datafile[, groups[i]] == as.character(extraction[[i]][j, 1])], na.rm = T), decimals) else extraction[[i]][j, proportion] <- round(mean(datafile[, response][datafile[, split.term(groups[i])[1]] == split.term(as.character(extraction[[i]][j, 1]))[1] & datafile[, split.term(groups[i])[2]] == split.term(as.character(extraction[[i]][j, 1]))[2]]), decimals)} if (fam == "binomial" & gold){ weight <- round(get.weights(extraction[[i]]$coefs, extraction[[i]]$tokens), decimals) weight[weight == 1] <- paste(" > 0.", paste(rep(9, decimals), collapse = ""), sep = "") weight[weight == 0] <- paste(" < 0.", paste(rep(0, decimals-1),collapse=""),1,sep = "") extraction[[i]]$weight <- as.character(weight) if (centered) names(extraction[[i]])[which(names(extraction[[i]]) == "weight")] <- "centered factor weight" else names(extraction[[i]])[which(names(extraction[[i]]) == "weight")] <- "uncentered factor weight"}} extraction[[i]]$coefs <- round(extraction[[i]]$coefs, decimals) if (length(extraction[[i]][, 1]) > 1) extraction[[i]] <- extraction[[i]][order(extraction[[i]]$coefs, decreasing = T), ] if (fam == "binomial") names(extraction[[i]])[2] <- "logodds" if (fam == "gaussian") names(extraction[[i]])[2] <- "coef" }} v <- try(vif.dej(model), silent = T) if (!"try-error" %in% class(v)){ if (all(is.numeric(v))) { for (x in names(extraction)){ if (v[x] > 2.5) f <- 2.5 else next if (v[x] > 5) f <- 5 if (v[x] > 7.5) f <- 7.5 if (v[x] > 10) f <- 10 if (v[x] > 15) f <- 15 if (v[x] > 20) f <- 20 extraction[[x]]$vif[1] <- paste(">", f, sep = "")}}} repeat { if ("glm" %in% class(model)) break r.den <- names(VarCorr(model)) for (r.d in r.den){ std.dev <- attr(VarCorr(model)[[r.d]], "stddev") names(std.dev)[1] <- "intercept" std.dev <- round(std.dev, decimals) extraction[[length(extraction) + 1]] <- as.data.frame(cbind(" " = "std dev", rbind(std.dev)), stringsAsFactors = F) names(extraction)[[length(extraction)]] <- paste(r.d, "(random)") ran <- as.character(eval(parse(text = paste(factor.slot, r.d, sep = "")))) extraction[[length(extraction)]]$tokens <- length(ran) if (fam == "gaussian") extraction[[length(extraction)]][proportion] <- round(mean(datafile[, response]), decimals) if (fam == "binomial") extraction[[length(extraction)]][proportion] <- round(mean(datafile[, response] == apps), decimals) } if (blups == F) break for (ran.d in paste(r.den, "(random)")){ r.d <- substr(ran.d, 1, nchar(ran.d) - 9) extraction[[ran.d]] <- rbind(extraction[[ran.d]], "...") namess <- levels(eval(parse(text=paste(factor.slot, r.d, sep = "")))) ran <- as.character(eval(parse(text = paste(factor.slot, r.d, sep = "")))) ranefs <- cbind(namess, round(as.data.frame(ranef(model)[r.d]), decimals), xtabs(~ ran)[namess]) if (fam == "gaussian") ranefs <- cbind(ranefs, as.vector(round(xtabs(eval(parse(text = response)) ~ eval(parse(text = r.d)), datafile) / xtabs(~ ran), decimals))) if (fam == "binomial") ranefs <- cbind(ranefs, as.vector(round(xtabs(~ eval(parse(text = r.d)), datafile[datafile[, response] == apps, ]) / xtabs(~ eval(parse(text = r.d)), datafile), decimals))) names(ranefs) <- names(extraction[[ran.d]]) if (fam == "binomial" & gold & length(random.s[grep(r.d, random.s)]) == 0){ if (centered) fw <- "centered factor weight" else fw <- "uncentered factor weight" ranefs[, fw] <- round(get.weights(ranefs[, 2], xtabs(~ ran)[namess]), decimals) extraction[[ran.d]][, fw] <- "..."} extraction[[ran.d]] <- rbind(extraction[[ran.d]], ranefs[order(ranefs$intercept, decreasing = T), ]) } break } repeat { if ("glm" %in% class(model)) break if (contr == "contr.treatment") break r.den <- names(VarCorr(model)) for (r.d in r.den){ ran.d <- paste(r.d, "(random)") g <- random.s[grep(r.d, random.s)] slopes <- substr(g, 1, nchar(g) - nchar(r.d) - 1) for (s in slopes){ if (class(eval(parse(text=paste(factor.slot, s, sep = "")))) != "factor") next levs <- levels(eval(parse(text=paste(factor.slot, s, sep = "")))) oldcols <- names(extraction[[ran.d]]) minbro <- min(grep(s, oldcols)) maxbro <- max(grep(s, oldcols)) zs <- "..." L <- length(extraction[[ran.d]][, 1]) if (L > 2){ zs <- c("NA", "...") for (i in 3:L){ zs <- c(zs, round(-sum(as.numeric( extraction[[ran.d]][i, minbro:maxbro])), decimals))}} gg <- grep(s, colnames(VarCorr(model)[[r.d]])) vc <- VarCorr(model)[[r.d]][gg, gg] zs[1] <- round(sqrt(sum(vc)), decimals) newsection <- cbind(extraction[[ran.d]][, 1:maxbro], data.frame(zs = zs)) if (maxbro < length(oldcols)) newsection <- cbind(newsection, extraction[[ran.d]][, (maxbro + 1):length(oldcols)]) extraction[[ran.d]] <- newsection oldcols[minbro:maxbro] <- paste(s, levs[1:(length(levs) - 1)], sep = "") newcol <- paste(s, levs[length(levs)], sep = "") names(extraction[[ran.d]]) <- c(oldcols[1:maxbro], newcol) if (maxbro < length(oldcols)) names(extraction[[ran.d]]) <- c(oldcols[1:maxbro], newcol, oldcols[(maxbro + 1):length(oldcols)])}} break } extraction$misc.1 <- data.frame(n = n, df = "", intercept = intercept, mean = "") if (fam == "binomial" & gold & all(!names(extraction) %in% continuous) & !all(names(extraction) %in% c("misc.1", "misc.2"))){ if (centered) extraction$misc.1$"centered input prob" <- round(get.input(extraction), decimals) else extraction$misc.1$"uncentered input prob" <- round(get.input(extraction), decimals)} if ("glm" %in% class(model)){ slashes <- length(which(abs(coef(model)) > 0 & abs(coef(model)) < 1e-50)) extraction$misc.1$df <- length(datafile[, response]) - model$df.residual - slashes} if (class(model)[1] %in% c("lmerMod", "glmerMod")){ slashes <- length(which(abs(fixef(model)) > 0 & abs(fixef(model)) < 1e-50)) extraction$misc.1$df <- attr(logLik(model), "df") - slashes} extraction$misc.1$intercept <- round(intercept, decimals) if (fam == "binomial"){ extraction$misc.1$mean <- round(length(datafile[datafile[, response] == apps, response])/ (length(datafile[datafile[, response] == apps, response]) + length(datafile[datafile[, response] == non.apps, response])), decimals) names(extraction$misc.1)[which(names(extraction$misc.1) == "mean")] <- "overall proportion"} if (fam == "gaussian"){ extraction$misc.1$mean <- round(mean(datafile[, response], na.rm = T), decimals) names(extraction$misc.1)[which(names(extraction$misc.1) == "mean")] <- "overall mean"} extraction$misc.2 <- data.frame(deviance = "") if ("glm" %in% class(model)){ extraction$misc.2$deviance <- round(model$deviance, decimals) extraction$misc.2$AIC <- round(AIC(model), decimals) extraction$misc.2$AICc <- round(AICc(model), decimals) if (fam == "binomial"){ extraction$misc.2$Dxy <- round(unname(somers2(fitted(model), model$y)[2]), decimals)} r2 <- try(round(suppressMessages(r.squaredGLMM(model)), decimals), silent = T) if (!"try-error" %in% class(r2)){ extraction$misc.2$R2 <- r2[2]}} if (class(model)[1] %in% c("lmerMod", "glmerMod")){ extraction$misc.2$deviance <- round(deviance(model, REML = F), decimals) extraction$misc.2$AIC <- round(AIC(model), decimals) extraction$misc.2$AICc <- round(AICc(model), decimals) if (fam == "binomial"){ extraction$misc.2$Dxy.fixed <- unname(round(somers2(plogis(predict(model, re.form = NA)), model@resp$y)[2]), decimals) extraction$misc.2$Dxy.total <- round(unname(somers2(fitted(model), model@resp$y)[2]), decimals)} options(warn = -1) r2 <- try(round(suppressMessages(r.squaredGLMM(model)), decimals), silent = T) options(warn = 2) if (!"try-error" %in% class(r2)){ extraction$misc.2$R2.fixed <- r2[1] extraction$misc.2$R2.random <- r2[2] - r2[1] extraction$misc.2$R2.total <- r2[2]} } extraction} fit.model <- function(data = datafile){ if (is.null(contr)) contr <<- "contr.sum" options(contrasts = c(unordered = contr, ordered = "contr.poly")) if (length(fixed) == 0) fixed <<- 1 if (length(ff <- fixed[fixed != 1]) > 0){ old.contrasts <- list() for (f in ff){ num <- which(ff == f) old.contrasts[[num]] <- contrasts(datafile[, f]) names(old.contrasts)[num] <- f if (!"/" %in% levels(datafile[, f])) next if (length(levels(datafile[, f])) < 3) next firstlev <- levels(datafile[, f])[1] if (firstlev == "/") firstlev <- levels(datafile[, f])[2] datafile[, f] <<- relevel(datafile[, f], "/") if (contr == "contr.treatment") datafile[, f] <<- relevel(datafile[, f], firstlev) if (contr == "contr.sum") contrasts(datafile[, f])[1, 1] <<- 1e+100 if (contr == "contr.treatment") contrasts(datafile[, f])[2, 1] <<- 1e+100}} if (length(continuous) == 0) continuous <<- 1 if (length(interactions) == 0) interactions <<- 1 form <- paste(response, paste( paste(fixed, collapse = "+"), paste(continuous, collapse = "+"), paste(interactions, collapse = "+"), sep = "+"), sep = "~") if (length(random.i) > 0){ for (r in random.i){ slope.full <- random.s[grep(paste(".", r, sep = ""), random.s)] slope <- paste(substr(slope.full, 1, nchar(slope.full) - nchar(r) - 1), collapse = "+") if (slope == "") slope <- 1 form <- paste(form, " + (", slope, "|", r, ")", sep = "")} if (fam == "gaussian"){ model <- try(lmer(as.formula(form), data, control = lmerControl(optCtrl = list(maxfun = maxfun))), silent = T) if ("try-error" %in% class(model)){ if (length(grep("failed to converge", model) > 0)){ options(warn = 1) model <- try(lmer(as.formula(form), data, control = lmerControl(check.conv.grad = .makeCC("warning", tol = tolerance, relTol = NULL), optCtrl = list(maxfun = maxfun))), silent = T) options(warn = 2)}}} if (fam == "binomial"){ model <- try(glmer(as.formula(form), data, fam, control = glmerControl(optCtrl = list(maxfun = maxfun))), silent = T) if ("try-error" %in% class(model)){ if (length(grep("failed to converge", model) > 0)){ options(warn = 1) model <- try(glmer(as.formula(form), data, fam, control = glmerControl(check.conv.grad = .makeCC("warning", tol = tolerance, relTol = NULL), optCtrl = list(maxfun = maxfun))), silent = T) options(warn = 2)}}}} else model <- try(glm(as.formula(form), fam, data), silent = T) for (f in ff) contrasts(datafile[, f]) <<- old.contrasts[[f]] model} get.contrasts <- function(x, n){ cont <- vector() for (i in 1:length(x)) cont[i] <- contrasts(x)[x[i], n] cont} get.continuous <- function(){ continuous <<- vector() cat(sep = "\n", paste("Are any predictors continuous? (", get.menu(predictors, predictors.pick, paren = F), " Enter-none)", sep = "")) repeat {continuous.pick <<- try(as.numeric(scan(what = "character", nmax = length(predictors), quiet = T)), silent = T) continuous <<- names(datafile)[continuous.pick] if (!NA %in% continuous & all(continuous %in% predictors)) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(continuous) > 0) for (i in 1:length(continuous)){ contin <- try(as.numeric(as.character(datafile[, continuous[i]])), silent = T) if (!"try-error" %in% class(contin)) datafile[, continuous[i]] <<- contin else cat(sep = "\n", paste("Warning: could not make", continuous[i], "continuous!"))} for (j in 1:length(predictors)) if (!predictors[j] %in% continuous) datafile[, predictors[j]] <<- as.factor(datafile[, predictors[j]])[drop = T]} get.data <- function(){ if (!is.null(datafile)) save.data() cat("\n") if (is.null(datafile)) cat(sep = "\n", "No data loaded.") else cat(sep = "\n", paste("Current data file is:", fullpath)) cat("\n") cat(sep = "\n", "What separates the columns in the data file to open?") cat(sep = "\n", "(c-commas s-semicolons t-tabs tf-token file)") cat(sep = "\n", "Press Enter to exit, keeping current data file, if any.") repeat{sep.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(sep.pick) == 0) return() if (sep.pick %in% c("c", "s", "t", "tf")) break} if (sep.pick == "c"){ datafile <<- try(read.table(fullpath <<- file.choose(), header = T, sep = ",", quote = "\"", comment.char = ""), silent = T) sepchar <<- ","} if (sep.pick == "s"){ datafile <<- try(read.table(fullpath <<- file.choose(), header = T, sep = ";", quote = "\"", comment.char = ""), silent = T) sepchar <<- ";"} if (sep.pick == "t"){ datafile <<- try(read.table(fullpath <<- file.choose(), header = T, sep = "\t", quote = "\"", comment.char = ""), silent = T) sepchar <<- "\t"} if (sep.pick == "tf"){ raw <- try(readLines(fullpath <<- file.choose()), silent = T) if ("try-error" %in% class(raw)) datafile <<- raw else { datlines <- grep("^\\(", raw) dat <- substr(raw[datlines], 2, nchar(raw[datlines])) spaces <- regexpr(" +", dat) comm <- substr(dat, spaces + attr(spaces, "match.length"), nchar(dat)) rump <- substr(dat, 1, ifelse(spaces > -1, spaces-1, nchar(dat))) minw <- min(nchar(rump)) maxw <- max(nchar(rump)) if (minw!= maxw){ cat(sep = "\n", "Warning! Tokens not all the same length. Omit short token(s) or drop final column(s)?") cat(sep = "\n", "(o-omit tokens d-drop columns)") repeat{drop.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(drop.pick) > 0) if (drop.pick %in% c("d", "o")) break} if (drop.pick == "o"){ rump <- rump[which(nchar(rump) == maxw)] minw <- maxw}} datafile <<- data.frame(dummy = 1:length(rump)) for (i in 1:minw) datafile[, i+1] <<- as.factor(substr(rump, i, i)) datafile <<- datafile[, -1] sepchar <<- ", " names(datafile) <<- 1:minw cat("\n") print(datafile[1:5, ]) cat("\n") repeat {cat(sep = "\n", "Enter labels for the columns (press Enter twice when done)") names <- try(scan(what = "character", nmax = minw, quiet = T)) if (!"try-error" %in% names & length(names) > 0) break cat(sep = "\n", "Enter at least two legitimate column names!")} names(datafile) <<- names datafile <<- datafile[, 1:length(names(datafile)[!is.na(names(datafile))])] if (max(spaces) > -1) datafile$comments <<- ifelse(spaces > -1, comm, "")}} if ("try-error" %in% class(datafile)){ cat("\n") cat(sep = "\n", paste("Error in loading file.")) cat(sep = "\n", datafile[1]) rm(datafile, envir = globalenv()) return()} lastslash <- gregexpr("/", fullpath)[[1]][length(gregexpr("/", fullpath)[[1]])] if (lastslash == -1) lastslash <- gregexpr("\\\\", fullpath)[[1]][length(gregexpr("\\\\", fullpath)[[1]])] try(setwd(shortpath <<- substr(fullpath, 1, lastslash)), silent = T) droprow <- vector() for (i in 1:length(datafile[, 1])) if (all(is.na(datafile[i, ]))) droprow <- c(droprow, i) if (length(droprow) > 0) datafile <<- datafile[-droprow, ] dropcolumn <- vector() for (j in 1:length(datafile[1, ])) if (all(is.na(datafile[, j]))) dropcolumn <- c(dropcolumn, j) if (length(dropcolumn) > 0) datafile <<- datafile[, -dropcolumn] loaded.flag <<- T datafile2 <<- datafile } get.input <- function(extraction){ tally <- 0 if (length(extraction) > 1){ for (i in 1:(length(extraction) - 1)){ if (names(extraction[[i]])[1] == "factor") { e <- extraction[[i]] e$logodds[e$factor == "/"] <- 0 if (centered) tally <- tally + sum(e$logodds) / length(e$logodds) else tally <- tally + sum(e$logodds * e$tokens) / sum(e$tokens)}}} plogis(extraction$misc.1$intercept + tally)} get.interactions <- function(){ interactions <<- vector() interactions.picks <<- vector() interactions.update <- vector() if (length(predictors) == 1) return() for (i in 1:sum(1:length(predictors)-1)){ cat(sep = "\n", paste("Consider the interaction between two fixed effects? [For interactions between fixed and random effects, use random slopes, below.] (", get.menu(predictors, predictors.pick, paren = F), " Enter-done)", sep = "")) repeat {interactions.pick <- try(as.numeric(scan(what = "character", nmax = 2, quiet = T)), silent = T) interactions.update <- paste(names(datafile)[interactions.pick], collapse = ":") if (length(interactions.pick) == 2 & !"NA" %in% split.term(interactions.update) & !interactions.update %in% interactions) if (!identical(split.term(interactions.update)[1], split.term(interactions.update)[2])) break if (length(interactions.pick) == 0) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(interactions.pick) == 0) break interactions <<- c(interactions, interactions.update) interactions.picks <<- c(interactions.picks, interactions.pick)}} get.menu <- function(variable, index = 1:length(variable), paren = T){ if (paren) menu <- paste(" (", paste(index, "-", variable, sep = "", collapse = " "), ")", sep = "") if (!paren) menu <- paste(index, "-", variable, sep = "", collapse = " ") menu} get.p <- function(model1, model2){ if ("try-error" %in% class(model1) & "try-error" %in% class(model2)){ p.flag <<- "error code = " p <- ifelse(model2 == "interaction", -888, -999) return(p)} if ("try-error" %in% class(model1) & !"try-error" %in% class(model2)){ mod2 <- extract(model2) p.flag <<- "log-likelihood = " p <- -mod2$misc.2$deviance/2 return(p)} if (!"try-error" %in% class(model1) & "try-error" %in% class(model2)){ p.flag <<- "error code = " p <- ifelse(model2 == "interaction", 888, 999) return(p)} p.flag <<- "p = " if ("glm" %in% class(model1) & "glm" %in% class(model2)){ if (!identical(model1$model[, 1], model2$model[, 1])) return("error") deviance1 <- -2 * logLik(model1) deviance2 <- -2 * logLik(model2) df1 <- attr(deviance1, "df") - length(which(abs(coef(model1)) > 0 & abs(coef(model1)) < 1e-50)) df2 <- attr(deviance2, "df") - length(which(abs(coef(model2)) > 0 & abs(coef(model2)) < 1e-50)) if ((dd <- abs(deviance1 - deviance2)) == 0) return(1) p <- as.vector(pchisq(dd, abs(df1 - df2), lower = F)) if (is.na(p)) return(1) return(p)} if (class(model1)[1] %in% c("lmerMod", "glmerMod") & class(model2)[1] %in% c("lmerMod", "glmerMod")){ if (!identical(model1@resp$y, model2@resp$y)) return("error") options(warn = 1) model11 <- refitML(model1) model22 <- refitML(model2) deviance11 <- -2 * logLik(model11) deviance22 <- -2 * logLik(model22) df11 <- attr(deviance11, "df") - length(which(abs(fixef(model11)) > 0 & abs(fixef(model11)) < 1e-50)) df22 <- attr(deviance22, "df") - length(which(abs(fixef(model22)) > 0 & abs(fixef(model22)) < 1e-50)) options(warn = 2) if ((dd <- abs(deviance11 - deviance22)) == 0) return(1) p <- as.vector(pchisq(dd, abs(df11 - df22), lower = F)) if (is.na(p)) return(1) return(p)} "error"} get.packages <- function(packages = T){ repo <- options("repos") options(repos = "http://cran.r-project.org") if (packages == F){ try1a <- try(suppressPackageStartupMessages(library(lme4)), silent = T) try2a <- try(suppressPackageStartupMessages(library(MuMIn)), silent = T) try3a <- try(suppressPackageStartupMessages(library(Hmisc)), silent = T) try4a <- try(suppressPackageStartupMessages(library(car)), silent = T) try5a <- try(suppressPackageStartupMessages(library(lattice)), silent = T) return()} version <- paste(sessionInfo()$R.version$major, sessionInfo()$R.version$minor, sep = ".") if ("lme4" %in% installed.packages()[, 1]){ if (version < "3.2.2") options(warn = -1) try1 <- try(update.packages(oldPkgs = "lme4", ask = F, type = "binary"), silent = T) if (version < "3.2.2") options(warn = 2) if ("try-error" %in% class(try1)){ cat("\n") cat(paste("Problem updating lme4 package but package is installed:\n", try1, sep = ""))}} else { if (version < "3.2.2") options(warn = -1) try11 <- try(install.packages("lme4", type = "binary", quiet = T), silent = T) if (version < "3.2.2") options(warn = 2) if ("try-error" %in% class(try11)){ cat("\n") cat(paste("Problem installing lme4 package:\n", try11, sep = ""))}} try111 <- try(suppressPackageStartupMessages(library(lme4)), silent = T) if ("try-error" %in% class(try111)){ cat("\n") cat(paste("Problem loading lme4 package:\n", try111, sep = ""))} else { cat("\n") cat("Successfully loaded lme4 package.")} if ("MuMIn" %in% installed.packages()[, 1]){ try2 <- try(update.packages(oldPkgs = "MuMIn", ask = F, type = "binary"), silent = T) if ("try-error" %in% class(try2)){ cat("\n") cat(paste("Problem updating MuMIn package but package is installed:\n", try2, sep = ""))}} else { try22 <- try(install.packages("MuMIn", type = "binary", quiet = T), silent = T) if ("try-error" %in% class(try22)){ cat("\n") cat(paste("Problem installing MuMIn package:\n,", try22, sep = ""))}} try222 <- try(suppressPackageStartupMessages(library(MuMIn)), silent = T) if ("try-error" %in% class(try222)){ cat("\n") cat(paste("Problem loading MuMIn package:\n", try222, sep = ""))} else { cat("\n") cat("Successfully loaded MuMIn package.")} if ("Hmisc" %in% installed.packages()[, 1]){ try3 <- try(update.packages(oldPkgs = "Hmisc", ask = F, type = "binary"), silent = T) if ("try-error" %in% class(try3)){ cat("\n") cat(paste("Problem updating Hmisc package but package is installed:\n", try3, sep = ""))}} else { try33 <- try(install.packages("Hmisc", type = "binary", quiet = T), silent = T) if ("try-error" %in% class(try33)){ cat("\n") cat(paste("Problem installing Hmisc package:\n,", try33, sep = ""))}} try333 <- try(suppressPackageStartupMessages(library(Hmisc)), silent = T) if ("try-error" %in% class(try333)){ cat("\n") cat(paste("Problem loading Hmisc package:\n", try333, sep = ""))} else { cat("\n") cat("Successfully loaded Hmisc package.")} if ("car" %in% installed.packages()[, 1]){ try4 <- try(update.packages(oldPkgs = "car", ask = F, type = "binary"), silent = T) if ("try-error" %in% class(try4)){ cat("\n") cat(paste("Problem updating car package but package is installed:\n", try4, sep = ""))}} else { try44 <- try(install.packages("car", type = "binary", quiet = T), silent = T) if ("try-error" %in% class(try44)){ cat("\n") cat(paste("Problem installing car package:\n,", try44, sep = ""))}} try444 <- try(suppressPackageStartupMessages(library(car)), silent = T) if ("try-error" %in% class(try444)){ cat("\n") cat(paste("Problem loading car package:\n", try444, sep = ""))} else { cat("\n") cat("Successfully loaded car package.")} if ("lattice" %in% installed.packages()[, 1]){ try5 <- try(update.packages(oldPkgs = "lattice", ask = F, type = "binary"), silent = T) if ("try-error" %in% class(try5)){ cat("\n") cat(paste("Problem updating lattice package but package is installed:\n", try5, sep = ""))}} else { try55 <- try(install.packages("lattice", type = "binary", quiet = T), silent = T) if ("try-error" %in% class(try55)){ cat("\n") cat(paste("Problem installing lattice package:\n,", try55, sep = ""))}} try555 <- try(suppressPackageStartupMessages(library(lattice)), silent = T) if ("try-error" %in% class(try555)){ cat("\n") cat(paste("Problem loading lattice package:\n", try555, sep = ""))} else { cat("\n") cat("Successfully loaded lattice package.")} options(repos = repo)} get.predictors <- function(){ if (!is.null(predictors) & !loaded.flag) old.predictors <- predictors singletons <- vector() for (i in names(datafile)[-response.pick]) if (length(levels(datafile[, i])) == 1) singletons <- c(singletons, i) cat(sep = "\n", paste("Choose predictors (independent variables) by number", if (!is.null(old.predictors)) if (!response %in% old.predictors) paste(", or Enter to keep ", paste(old.predictors, collapse = " & "), sep = ""), get.menu(names(datafile)[!names(datafile) %in% c(response, singletons)], which(!names(datafile) %in% c(response, singletons))), sep = "")) repeat {predictors.pick <<- try(as.numeric(scan(what = "character", nmax = length(names(datafile))-length(c(response, singletons)), quiet = T)), silent = T) if (length(predictors.pick) == 0 & !is.null(old.predictors)){ if (!response %in% old.predictors){ predictors <<- old.predictors predictors.pick <<- which(names(datafile) %in% predictors) return(T)}} predictors <<- names(datafile)[predictors.pick] if (!NA %in% predictors & all(!predictors %in% c(response, singletons)) & length(predictors) > 0) break if (length(predictors) == 0) cat(sep = "\n", "Please choose at least one predictor.") else cat(sep = "\n", "Invalid entry, please try again.")} for (i in 1:length(predictors)) if (NA %in% datafile[, predictors[i]]){ cat("\n") cat(sep = "\n", paste("Warning! NA's exist in ", predictors[i], "! This may cause a crash if you try use ", predictors[i], "!", sep = "")) cat("\n")} if ("/" %in% as.character(datafile[, predictors[i]])){ cat("\n") cat(sep = "\n", paste("Slashes -> / <- exist in ", predictors[i], ". Use to 'ignore' tokens in calculating that predictor's effects.", sep = "")) cat("\n")} loaded.flag <<- F F} get.random <- function(){ random.i <<- vector() categorical <- predictors[!predictors %in% continuous] if (length(interactions) > 0){ interactants <- vector() for (int in interactions) interactants <- c(interactants, split.term(int)) categorical <- categorical[!categorical %in% interactants]} if (length(categorical) > 0){ cat(sep = "\n", paste("Any random intercepts? (", get.menu(categorical, predictors.pick[!predictors.pick %in% c(continuous.pick, interactions.picks)], paren = F), " Enter-none)", sep = "")) repeat {random.i.pick <- try(as.numeric(scan(what = "character", nmax = length(categorical), quiet = T)), silent = T) random.i <<- names(datafile)[random.i.pick] if (!NA %in% random.i & all(random.i %in% categorical)) break cat(sep = "\n", "Invalid entry, please try again.")}} random.s <<- vector() potential <- c(predictors, interactions) potential <- potential[!potential %in% random.i] if (length(potential) > 0){ for (r in random.i){ cat(sep = "\n", paste("Any by-", r, " random slopes - must vary for each ", r, "? (", get.menu(potential[potential != r], which(names(datafile) %in% potential[potential != r]), paren = F), " Enter-none)", sep = "")) repeat {random.s.pick <- try(as.numeric(scan(what = "character", nmax = length(potential), quiet = T)), silent = T) random.s.update <- vector() if (length(random.s.pick) == 0) break random.s.update <- paste(names(datafile)[random.s.pick], r, sep = ".") if (!NA %in% random.s.update & all(random.s.update %in% paste(potential, r, sep = "."))) break cat(sep = "\n", "Invalid entry, please try again.")} random.s <<- c(random.s, random.s.update)}}} get.response <- function(){ if (!is.null(response) & !loaded.flag) old.response <- response singletons <- vector() names(datafile)[which(names(datafile) == "response")] <<- "response." for (i in names(datafile)) if (length(levels(datafile[, i])) == 1) singletons <- c(singletons, i) cat(sep = "\n", paste("Choose response (dependent variable) by number", if (!is.null(old.response)) paste(", or Enter to keep ", old.response, sep = ""), get.menu(names(datafile)[!names(datafile) %in% singletons], which(!names(datafile) %in% singletons)), sep = "")) repeat {response.pick <<- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(response.pick) == 0 & !is.null(old.response)) { response <<- old.response response.pick <<- which(names(datafile) == response) return(list(response, datafile[, response]))} response <<- names(datafile)[response.pick] if (!NA %in% response & length(response) > 0) break if (length(response) == 0) cat(sep = "\n", "Please choose a response variable.") else cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Type of response? (1-continuous Enter-binary)") repeat {fam.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(fam.pick) == 0) {fam <<- "binomial"; break} if (fam.pick == 1) {fam <<- "gaussian"; break} cat(sep = "\n", "Invalid entry, please try again.")} if (fam == "binomial"){ datafile[, response] <<- as.factor(datafile[, response])[drop = T] fac <- levels(datafile[, response]) cat(sep = "\n", paste("Choose application value(s) by number?", get.menu(fac))) repeat {application.values.pick <- try(as.numeric(scan(what = "character", nmax = length(fac)-1, quiet = T)), silent = T) application.values <- fac[application.values.pick] if (length(application.values) > 0) if (!NA %in% application.values) break cat(sep = "\n", "Invalid entry, please try again.")} apps <<- paste(application.values, collapse = "+") for (i in 1:length(fac)) if (fac[i] %in% application.values) fac[i] <- apps if (length(fac[fac != apps]) == 1) non.application.values <- fac[fac != apps] else { cat(sep = "\n", paste("Choose non-application value(s) by number, or press Enter for all remaining factors?")) cat(sep = "\n", get.menu(fac[fac!= apps], which(fac!= apps))) repeat {non.application.values.pick <- try(as.numeric(scan(what = "character", nmax = length(fac)-length(application.values), quiet = T)), silent = T) non.application.values <- fac[non.application.values.pick] if (!NA %in% non.application.values & all(non.application.values %in% fac[fac!= apps])) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(non.application.values) == 0) non.application.values <- fac[fac != apps]} non.apps <<- paste(non.application.values, collapse = "+") proportion <<- paste(apps, "/", apps, "+", non.apps, sep = "") for (i in 1:length(fac[fac != apps])) if (fac[fac != apps][i] %in% non.application.values) fac[fac!= apps][i] <- non.apps levels(datafile[, response]) <<- fac datafile[, response] <<- relevel(datafile[, response], non.apps) datafile <<- subset(datafile, datafile[, response] %in% c(non.apps, apps, NA)) datafile[, response] <<- datafile[, response][drop = T]} if (fam == "gaussian"){ proportion <<- "mean" resp <- try(as.numeric(as.character(datafile[, response])), silent = T) if (!"try-error" %in% class(resp)) datafile[, response] <<- resp else {cat(sep = "\n", paste("Could not make", response, "continuous!")); get.response()}} if (NA %in% datafile[, response]){ cat("\n") cat(sep = "\n", "Warning! NA's exist in response! This may not cause a crash, but it may not be what you want! Observations where response is NA will be ignored in all models.") cat("\n")} list(response, datafile[, response])} get.settings <- function(){ cat(sep = "\n", "Run silently? (1-yes Enter-no)") repeat {silent.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(silent.pick) == 0) {silent <<- F; break} if (silent.pick == 1) {silent <<- T; break} cat(sep = "\n", "Invalid entry, please try again.")} if (!silent){ cat(sep = "\n", "Hide intermediate models' details? (1-yes Enter-no)") repeat {verbose.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(verbose.pick) == 0) {verbose <<- T; break} if (verbose.pick == 1) {verbose <<- F; break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Hide factor weights? (1-yes Enter-no)") repeat {gold.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(gold.pick) == 0) {gold <<- T; break} if (gold.pick == 1) {gold <<- F; centered <<- F; break} cat(sep = "\n", "Invalid entry, please try again.")} if (gold){ cat(sep = "\n", "Center factor weights [recommended]? (1-no Enter-yes)") repeat {centered.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(centered.pick) == 0) {centered <<- T; break} if (centered.pick == 1) {centered <<- F; break} cat(sep = "\n", "Invalid entry, please try again.")}} cat(sep = "\n", "Show random effect estimates (BLUPs)? (1-no Enter-yes)") repeat {blups.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(blups.pick) == 0) {blups <<- T; break} if (blups.pick == 1) {blups <<- F; break} cat(sep = "\n", "Invalid entry, please try again.")}} cat(sep = "\n", "Threshold p-value for fixed effect significance? (1-automatic, Enter-0.05, or type another value)") repeat {thresh.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(thresh.pick) == 0) {thresh <<- 0.05; break} if (0 < thresh.pick & thresh.pick < 1) {thresh <<- thresh.pick; break} if (thresh.pick == 1) {thresh <<- "A"; break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Number of significant digits/decimal places to display? (Enter-3)") repeat {decimals.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(decimals.pick) == 0) {decimals <<- 3; break} if (!"try-error" %in% class(decimals.pick)) {decimals <<- decimals.pick; break}} cat(sep = "\n", "Type of contrasts to use for factors? (1-treatment Enter-sum)") repeat {contr.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(contr.pick) == 0) {contr <<- "contr.sum"; break} if (contr.pick == 1) {contr <<- "contr.treatment"; break} cat(sep = "\n", "Invalid entry, please try again.")} options(contrasts = c(unordered = contr, ordered = "contr.poly")) cat(sep = "\n", "Maximum number of function evaluations for mixed-effects models? (Enter-100000)") repeat {maxfun.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(maxfun.pick) == 0) {maxfun <<- 100000; break} if (class(maxfun.pick) != "try-error") {maxfun <<- maxfun.pick; break} cat(sep = "\n", "Invalid entry, please try again.")} cat(sep = "\n", "Tolerance for the gradient for model convergence (higher is less strict)? (Enter-0.002)") repeat {tolerance.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(tolerance.pick) == 0) {tolerance <<- 0.002; break} if (class(tolerance.pick) != "try-error") {tolerance <<- tolerance.pick; break} cat(sep = "\n", "Invalid entry, please try again.")}} get.variables <- function(){ r <- get.response() response <<- r[[1]] datafile[, response] <<- r[[2]] g <- get.predictors() if (!g){ get.continuous() get.interactions() get.random()} set.variables() load.variables() variables} get.weights <- function(logodds, tokens = rep(1, length(logodds))){ if (length(which(!is.na(logodds))) < 2) weights <- NA else { real.logodds <- logodds[which(!is.na(logodds))] if (real.logodds[1] != 0) real.logodds <- real.logodds-real.logodds[1] real.tokens <- tokens[which(!is.na(logodds))] num.weights <- length(real.logodds) right <- real.logodds if (centered) left <- rep(1/num.weights, num.weights) else left <- real.tokens/num.weights for (i in 1:(num.weights-1)){ left <- rbind(left, c(-1, rep(0, i-1), 1, rep(0, num.weights-1-i)))} real.weights <- plogis(solve(as.matrix(left), right)) weights <- rep(NA, length(logodds)) weights[which(!is.na(logodds))] <- real.weights} weights} load.variables <- function(){ response <<- variables[[1]] if (names(variables)[1] == "response.binary") fam <<- "binomial" if (names(variables)[1] == "response.continuous") fam <<- "gaussian" fixed <<- vector() continuous <<- vector() interactions <<- vector() random.i <<- vector() random.s <<- vector() for (i in 2:length(variables)){ if (names(variables)[i] == "fixed.factor") fixed <<- c(fixed, variables[[i]]) if (names(variables)[i] == "fixed.continuous") continuous <<- c(continuous, variables[[i]]) if (names(variables)[i] == "fixed.interaction") interactions <<- c(interactions, variables[[i]]) if (names(variables)[i] == "random.intercept") random.i <<- c(random.i, variables[[i]]) if (names(variables)[i] == "random.slope") random.s <<- c(random.s, variables[[i]])}} modeling <- function(){ modeling.pick <<- 99 while (!modeling.pick %in% c(9, 0)){ cat("\n") if (is.null(variables)) {cat(sep = "\n", "No variables chosen."); cat("\n")} else { cat(sep = "\n", "Current variables are:") for (i in 1:length(variables)){ cat(paste(names(variables)[i], ": ", paste(variables[[i]], collapse = " "), sep = "")) if (names(variables)[i] == "response.binary") cat(paste(" (", apps, " vs. ", non.apps, ")", sep = "")) cat("\n")} cat("\n")} cat(sep = "\n", "MODELING MENU") cat(sep = "\n", "1-choose variables 2-one-level (recommended)") cat(sep = "\n", "3-step-up 4-step-down 5-step-up/step-down") cat(sep = "\n", "6-trim 7-plotting 8-settings 9-main menu 0-exit") cat(sep = "\n", "10-chi-square test") modeling.pick <<- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(modeling.pick) == 0) modeling() if (modeling.pick == 1) get.variables() if (modeling.pick == 2) step.one() if (modeling.pick == 3) step.up() if (modeling.pick == 4) step.down() if (modeling.pick == 5) step.updown() if (modeling.pick == 6) trim() if (modeling.pick == 7) plotting() if (modeling.pick == 8) get.settings() if (modeling.pick == 10) test.chisq() if (modeling.pick == 9) plotting.pick <<- 9 if (modeling.pick == 0) {plotting.pick <<- 0; main.pick <<- 0}}} overwrite <- function(x, y){ x <- as.character(x) if (!exists(x)) {assign(x, y, envir = globalenv()); return()} if (length(eval(parse(text = x))) == 0) {assign(x, y, envir = globalenv()); return()} if (is.null(x)) {assign(x, y, envir = globalenv()); return()}} panel.rsuperpose <- function(x, y = NULL, subscripts, groups, panel.groups = panel.rxyplot, col = NA, col.line = superpose.line$col, col.symbol = superpose.symbol$col, pch = superpose.symbol$pch, cex = superpose.symbol$cex, fill = superpose.symbol$fill, font = superpose.symbol$font, lty = superpose.line$lty, lwd = superpose.line$lwd, alpha = superpose.symbol$alpha, type = "p", ..., distribute.type = FALSE, scale.factor = 0) { if (distribute.type) { type <- as.list(type) } else { type <- unique(type) wg <- match("g", type, nomatch = NA) if (!is.na(wg)) { panel.grid(h = -1, v = -1) type <- type[-wg] } type <- list(type) } x <- as.numeric(x) if (!is.null(y)) y <- as.numeric(y) if (length(x) > 0) { if (!missing(col)) { if (missing(col.line)) col.line <- col if (missing(col.symbol)) col.symbol <- col } superpose.symbol <- trellis.par.get("superpose.symbol") superpose.line <- trellis.par.get("superpose.line") vals <- if (is.factor(groups)) levels(groups) else sort(unique(groups)) nvals <- length(vals) col <- rep(col, length = nvals) col.line <- rep(col.line, length = nvals) col.symbol <- rep(col.symbol, length = nvals) pch <- rep(pch, length = nvals) fill <- rep(fill, length = nvals) lty <- rep(lty, length = nvals) lwd <- rep(lwd, length = nvals) alpha <- rep(alpha, length = nvals) cex <- rep(cex, length = nvals) font <- rep(font, length = nvals) type <- rep(type, length = nvals) panel.groups <- if (is.function(panel.groups)) panel.groups else if (is.character(panel.groups)) get(panel.groups) else eval(panel.groups) subg <- groups[subscripts] ok <- !is.na(subg) for (i in seq_along(vals)) { id <- ok & (subg == vals[i]) typehack <- type[[i]] if (vals[i] == "XXX"){ typehack[which(typehack == "p")] <- "pm" typehack[which(typehack == "l")] <- "lm"} if (any(id)) { args <- list(x = x[id], subscripts = subscripts[id], pch = pch[i], cex = cex[i], font = font[i], col = col[i], col.line = col.line[i], col.symbol = col.symbol[i], fill = fill[i], lty = lty[i], lwd = lwd[i], alpha = alpha[i], type = typehack, group.number = i, scale.factor = scale.factor, ...) if (!is.null(y)) args$y <- y[id] do.call(panel.groups, args) } } } } panel.rxyplot <- function (x, y, type = "p", groups = NULL, pch = if (is.null(groups)) plot.symbol$pch else superpose.symbol$pch, col, col.line = if (is.null(groups)) plot.line$col else superpose.line$col, col.symbol = if (is.null(groups)) plot.symbol$col else superpose.symbol$col, font = if (is.null(groups)) plot.symbol$font else superpose.symbol$font, lty = if (is.null(groups)) plot.line$lty else superpose.line$lty, cex = if (is.null(groups)) plot.symbol$cex else superpose.symbol$cex, fill = if (is.null(groups)) plot.symbol$fill else superpose.symbol$fill, lwd = if (is.null(groups)) plot.line$lwd else superpose.line$lwd, horizontal = FALSE, ..., jitter.x = FALSE, jitter.y = FALSE, factor = 0.5, amount = NULL, scale.factor = 0) { if (all(is.na(x) | is.na(y))) return() x <- as.numeric(x) y <- as.numeric(y) xx <- x[!is.na(y)] yy <- y[!is.na(y)] y.avg <- tapply(yy, xx, mean) xvals <- try(as.numeric(names(y.avg)), silent = T) if ("try-error" %in% class(xvals)) xvals <- as.factor(names(y.avg)) ord <- order(xvals) ym <- y.avg[ord] xm <- xvals[ord] if (scale.factor == 0) {nn <- 1; s <- 1} else {nn <- tapply(yy, xx, length); s <- scale.factor} plot.symbol <- trellis.par.get("plot.symbol") plot.line <- trellis.par.get("plot.line") superpose.symbol <- trellis.par.get("superpose.symbol") superpose.line <- trellis.par.get("superpose.line") if (!missing(col)) { if (missing(col.line)) col.line <- col if (missing(col.symbol)) col.symbol <- col } if (!is.null(groups)) panel.rsuperpose(x, y, type = type, groups = groups, pch = pch, col.line = col.line, col.symbol = col.symbol, font = font, lty = lty, cex = cex, fill = fill, lwd = lwd, horizontal = horizontal, panel.groups = panel.rxyplot, jitter.x = jitter.x, jitter.y = jitter.y, factor = factor, amount = amount, scale.factor = scale.factor, ...) else { if ("o" %in% type || "b" %in% type) type <- c(type, "p", "l") if ("g" %in% type) panel.grid(h = -1, v = -1) if ("p" %in% type) panel.points(x = if (jitter.x) jitter(x, factor = factor, amount = amount) else x, y = if (jitter.y) jitter(y, factor = factor, amount = amount) else y, cex = cex, fill = fill, font = font, col = col.symbol, pch = pch, ...) if ("pm" %in% type) panel.points(x = if (jitter.x) jitter(xm, factor = factor, amount = amount) else xm, y = if (jitter.y) jitter(ym, factor = factor, amount = amount) else ym, cex = cex*sqrt(nn)*s, fill = fill, font = font, col = col.symbol, pch = pch, ...) if ("l" %in% type) panel.lines(x = x, y = y, lty = lty, col = col.line, lwd = lwd, ...) if ("lm" %in% type) panel.lines(x = xm, y = ym, lty = lty, col = col.line, lwd = lwd, ...) if ("h" %in% type) { if (horizontal) panel.lines(x = x, y = y, type = "H", lty = lty, col = col.line, lwd = lwd, ...) else panel.lines(x = x, y = y, type = "h", lty = lty, col = col.line, lwd = lwd, ...) } if ("s" %in% type) { ord <- if (horizontal) sort.list(y) else sort.list(x) n <- length(x) xx <- numeric(2 * n - 1) yy <- numeric(2 * n - 1) xx[2 * 1:n - 1] <- x[ord] yy[2 * 1:n - 1] <- y[ord] xx[2 * 1:(n - 1)] <- x[ord][-1] yy[2 * 1:(n - 1)] <- y[ord][-n] panel.lines(x = xx, y = yy, lty = lty, col = col.line, lwd = lwd, ...) } if ("S" %in% type) { ord <- if (horizontal) sort.list(y) else sort.list(x) n <- length(x) xx <- numeric(2 * n - 1) yy <- numeric(2 * n - 1) xx[2 * 1:n - 1] <- x[ord] yy[2 * 1:n - 1] <- y[ord] xx[2 * 1:(n - 1)] <- x[ord][-n] yy[2 * 1:(n - 1)] <- y[ord][-1] panel.lines(x = xx, y = yy, lty = lty, col = col.line, lwd = lwd, ...) } if ("r" %in% type) panel.lmline(x, y, col = col.line, lty = lty, lwd = lwd, ...) if ("smooth" %in% type) panel.loess(x, y, horizontal = horizontal, col = col.line, lty = lty, lwd = lwd, ...) if ("a" %in% type) panel.linejoin(x, y, horizontal = horizontal, lwd = lwd, lty = lty, col.line = col.line, ...) } } partial <- function(quiet = F, partialrees = c(NULL, NULL)){ cat(sep = "\n", paste("Choose two factors to make an interaction group, to use *along with* the original two?\n(press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {partial.pick <- try(as.numeric(scan(what = "character", nmax = 2, quiet = T)), silent = T) partialrees <- names(datafile)[partial.pick] if (length(partialrees) == 2 & !NA %in% partialrees) if (!identical(partialrees[1], partialrees[2])) break cat(sep = "\n", "Invalid entry, please try again.")} x <- datafile[, partialrees[1]] y <- datafile[, partialrees[2]] if (!(class(x) %in% c("factor", "character") & class(y) %in% c("factor", "character"))){ cat(sep = "\n", "Those are not both factors!") return()} fac.x <- levels(x) cat(sep = "\n", paste("Factor(s) from ", partialrees[1], " to cross with factor(s) from ", partialrees[2], "?", get.menu(fac.x), sep = "")) repeat {x.pick <- try(as.numeric(scan(what = "character", nmax = length(fac.x), quiet = T)), silent = T) x.cross <- fac.x[x.pick] if (!NA %in% x.cross) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(x.cross) == 0) return() fac.y <- levels(y) cat(sep = "\n", paste("Factor(s) from ", partialrees[2], " to cross with factor(s) from ", partialrees[1], "?", get.menu(fac.y), sep = "")) repeat {y.pick <- try(as.numeric(scan(what = "character", nmax = length(fac.y), quiet = T)), silent = T) y.cross <- fac.y[y.pick] if (!NA %in% y.cross) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(y.cross) == 0) return() cat(sep = "\n", paste("Name for this interaction, or press Enter to use", paste(paste(x.cross, collapse = ""), paste(y.cross, collapse = ""), sep = "_"), "?")) interaction.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(interaction.pick) == 0) interaction.name <- paste(paste(x.cross, collapse = ""), paste(y.cross, collapse = ""), sep = "_") else interaction.name <- interaction.pick cat(sep = "\n", paste("Name for the other levels, or press Enter to use 'other'?")) other.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(other.pick) == 0) other.name <- "other" else other.name <- other.pick datafile[, paste(partialrees, collapse = "_")] <<- ifelse(x %in% x.cross & y %in% y.cross, interaction.name, other.name)} plotting <- function(){ plotting.pick <<- 99 while (!plotting.pick %in% c(9, 0)){ cat("\n") cat(sep = "\n", "PLOTTING MENU") cat(sep = "\n", "1-custom scatterplot") if (!is.null(datafile)) cat(sep = "\n", "5-modeling 7-save plot 9-main menu 0-exit") else cat(sep = "\n", "7-save plot 9-main menu 0-exit") plotting.pick <<- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(plotting.pick) == 0) plotting.pick <<- 0 if (plotting.pick == 1) custom.plot() if (plotting.pick == 5) modeling() if (plotting.pick == 7) cat(sep = "\n", "Ha! Not yet implemented!") if (plotting.pick == 9) modeling.pick <<- 9 if (plotting.pick == 0) {modeling.pick <<- 0; main.pick <<- 0}}} reclass <- function(){ data.structure() cat(sep = "\n", paste("Change class of which variable?", get.menu(names(datafile)), sep = "")) repeat {reclassee.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) reclassee <- names(datafile)[reclassee.pick] if (!NA %in% reclassee) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(reclassee) == 0) return() repeat{ cat(sep = "\n", paste("Change", reclassee, "to which class? (f-factor c-continuous [integer/numeric])")) reclass.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(reclass.pick) == 0) return() if (reclass.pick == "f"|reclass.pick == "F"){ target <- try(as.factor(datafile[, reclassee]), silent = T) if (!"try-error" %in% class(target)) datafile[, reclassee] <<- target else cat(sep = "\n", "Failed to change class!") return()} if (reclass.pick == "c"|reclass.pick == "C"){ target <- try(as.numeric(as.character(datafile[, reclassee])), silent = T) if (!"try-error" %in% class(target)) datafile[, reclassee] <<- target else cat(sep = "\n", "Failed to change class!") return()} cat(sep = "\n", "Invalid entry, please try again.")}} recode <- function(){ cat(sep = "\n", paste("Factor group to recode? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {recode.group.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) recode.group <- names(datafile)[recode.group.pick] if (!NA %in% recode.group) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(recode.group) == 0) return() if (class(datafile[, recode.group])!= "factor"){ cat(sep = "\n", paste("Please change class of", recode.group, "to 'factor' before recoding.")) return()} fac <- levels(datafile[, recode.group]) repeat{ cat(sep = "\n", paste("Factor(s) of ", recode.group, " to recode together? (", get.menu(fac, paren = F), " Enter-done)", sep = "")) repeat {source.factors.pick <- try(as.numeric(scan(what = "character", nmax = length(fac), quiet = T)), silent = T) source.factors <- fac[source.factors.pick] if (!NA %in% source.factors) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(source.factors) == 0) break sf <- paste(source.factors, sep = "", collapse = " ") cat(sep = "\n", paste("Recode ", sf, " as what?", sep = "")) target.factor <- scan(what = "character", nmax = 1, quiet = T) if (length(target.factor) == 0) break for (i in 1:length(fac)) if (fac[i] %in% source.factors) fac[i] <- target.factor} cat(sep = "\n", "Recode to new column? (Yes-type new column name No-press Enter)") repeat {new.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(new.pick) == 0) {levels(datafile[, recode.group]) <<- fac; break} a <- try(datafile[, new.pick] <<- datafile[, recode.group], silent = T) if (!"try-error" %in% class(a)) {levels(datafile[, new.pick]) <<- fac; break} cat(sep = "\n", "Invalid column name.")}} relev <- function(){ cat(sep = "\n", paste("Factor group to relevel? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {relev.group.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) relev.group <- names(datafile)[relev.group.pick] if (!NA %in% relev.group) break if (class(datafile[, relev.group])!= "factor"){ cat(sep = "\n", paste("Please change class of", relev.group, "to 'factor' before releveling.")) break} cat(sep = "\n", "Invalid entry, please try again.")} if (length(relev.group) == 0) return() fac <- levels(datafile[, relev.group]) cat(sep = "\n", paste("For", relev.group, "current baseline is", fac[1])) cat(sep = "\n", paste("Choose new baseline (by number)?", get.menu(fac))) repeat {relev.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) baseline <- fac[relev.pick] if (!NA %in% baseline) break cat ("Invalid entry, please try again.")} if (length(baseline) == 0) return() datafile[, relev.group] <<- relevel(datafile[, relev.group], baseline)} rename <- function(){ cat(sep = "\n", paste("Variable(s) to rename? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {rename.group.pick <- try(as.numeric(scan(what = "character", nmax = length(names(datafile)), quiet = T)), silent = T) rename.group <- names(datafile)[rename.group.pick] if (!NA %in% rename.group) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(rename.group) == 0) return() for (group in rename.group){ cat(sep = "\n", paste("Type new name for ", group, "?", sep = "")) repeat{ new.name <- scan(what = "character", nmax = 1, quiet = T) t <- try((names(datafile)[which(names(datafile) == group)] <<- new.name), silent = T) if ("try-error" %in% class(t)) cat(sep = "\n", "Invalid entry, please try again.") else break}}} reset <- function(){ for (v in vars) try(rm(list = eval(v), envir = globalenv()), silent = T) source("http://www.danielezrajohnson.com/Rbrul.R") rb <- try(rbrul2(packages = packages)) if ("try-error" %in% class(rb)){ cat(sep = "\n", "\nSorry, there was an error and Rbrul crashed.") cat(sep = "\n", "Help is available by emailing danielezrajohnson@gmail.com.")} main.pick <<- 0} retain <- function(){ cat(sep = "\n", paste("Factor group to retain using?", get.menu(names(datafile)), sep = "")) repeat {retain.group.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) retain.group <- names(datafile)[retain.group.pick] if (!NA %in% retain.group) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(retain.group) == 0) return() if (class(datafile[, retain.group])!= "factor"){ cat(sep = "\n", paste("Please change class of", retain.group, "to 'factor' before retaining.")) return()} fac <- levels(datafile[, retain.group]) cat(sep = "\n", paste("Factors to retain from ", retain.group, "?", get.menu(fac), sep = "")) repeat {retained.factors.pick <- try(as.numeric(scan(what = "character", nmax = length(fac), quiet = T)), silent = T) retained.factors <- fac[retained.factors.pick] if (!NA %in% retained.factors) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(retained.factors) == 0) return() datafile <<- subset(datafile, datafile[, retain.group] %in% retained.factors) datafile[, retain.group] <<- datafile[, retain.group][drop = T]} rownames.off <- function (x, ..., digits = NULL, quote = FALSE, right = TRUE, row.names = FALSE) { n <- length(row.names(x)) if (length(x) == 0L) { cat(gettextf("data frame with 0 columns and %d rows\n", n)) } else if (n == 0L) { print.default(names(x), quote = FALSE) cat(gettext(" < 0 rows > (or 0-length row.names)\n")) } else { m <- as.matrix(format.data.frame(x, digits = digits, na.encode = FALSE)) if (!isTRUE(row.names)) dimnames(m)[[1L]] <- if (identical(row.names, FALSE)) rep.int("", n) else row.names print(m, ..., quote = quote, right = right) } invisible(x) } save.data <- function(){ cat("\n") if (is.null(datafile)) { cat(sep = "\n", "No data loaded.") cat("\n")} else cat(sep = "\n", paste("Current data file is:", fullpath)) data.structure() cat(sep = "\n", "Save current data to a file? (Yes-type filename No-press Enter)") repeat {save.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(save.pick) == 0) break a <- try(write.table(datafile, save.pick, row.names = F, sep = sepchar, quote = F), silent = T) if (!"try-error" %in% class(a)){ if (length(grep("/", save.pick)) > 0) fullpath <<- save.pick else fullpath <<- paste(shortpath, save.pick, sep = "") break} cat(sep = "\n", "Error in saving file.")}} set.variables <- function(){ type <- "continuous" if (length(fam) > 0) if (fam == "binomial") type <- "binary" variables <<- list(response = response, fixed.factor = predictors[!predictors %in% continuous & !predictors %in% random.i], fixed.continuous = continuous, fixed.interaction = interactions, random.intercept = random.i, random.slope = random.s) names(variables)[1] <<- paste("response", type, sep = ".") blanks <- vector() for (j in 1:length(variables)) if (length(variables[[j]]) == 0) blanks <- c(blanks, j) if (length(blanks) > 0) variables <<- variables[-blanks]} split.term <- function(interaction){ colon <- regexpr(":", interaction) result <- vector() index <- 1 result[index] <- substr(interaction, 1, colon - 1) if (nchar(result[1]) > 0) index <- 2 result[index] <- substr(interaction, colon + 1, nchar(interaction)) result} star <- function(){ cat(sep = "\n", paste("Choose two factors to make an interaction group, to use *instead* of the original two?\n(press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {star.pick <- try(as.numeric(scan(what = "character", nmax = 2, quiet = T)), silent = T) starrees <- names(datafile)[star.pick] if (length(starrees) == 2 & !NA %in% starrees) if (!identical(starrees[1], starrees[2])) break cat(sep = "\n", "Invalid entry, please try again.")} x <- datafile[, starrees[1]] y <- datafile[, starrees[2]] if (class(x) %in% c("factor", "character") & class(y) %in% c("factor", "character")) datafile[, paste(starrees, collapse = ".")] <<- as.factor(paste(x, y, sep = ".")) else cat(sep = "\n", "Those are not both factors!")} step.down <- function(){ if (is.null(variables)) return(cat(sep = "\n", "Please choose variables first.")) load.variables() fix <- fixed cont <- continuous int <- interactions in.preds <- c(fixed, continuous, interactions) keep.going <- T steppe <- 0 run <- 1 best.model <- fit.model() if ("try-error" %in% class(best.model)) error.count <- 1 else error.count <- 0 best.run <- 0 if (!silent){ cat("\n") cat(sep = "\n", "STEPPING DOWN...") cat("\n") cat(sep = "\n", paste("STEP 0 - Run 0 - full model with", paste(in.preds, collapse = " + "))) cat("\n") if (verbose) print(extract(best.model), row.names = F)} while (keep.going & length(in.preds) > 0){ steppe <- steppe + 1 if (thresh == "A") threshold <<- 0.05/length(in.preds) else threshold <<- thresh compare <- vector() mod <- list() mod.run <- vector() if (!silent) {cat("\n") cat(sep = "\n", paste("STEP", steppe, "- dropping from Run", best.run, "model with", paste(in.preds, collapse = " + ")))} if (length(in.preds) > 0) for (i in 1:length(in.preds)){ if (in.preds[i] %in% sapply(in.preds[in.preds %in% int], split.term)){ mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error" compare[i] <- get.p(best.model, mod[[i]]) next} cat("\n") cat(sep = "\n", paste("Trying without ", in.preds[i], "...", sep = "")) selection <- in.preds[-i] continuous <<- cont[cont %in% selection] fixed <<- fix[fix %in% selection] interactions <<- int[int %in% selection] mod[[i]] <- fit.model() mod.run[i] <- run compare[i] <- get.p(best.model, mod[[i]]) if (!silent) {cat("\n") if (verbose) print(extract(mod[[i]]), row.names = F) above.preds <- ifelse(length(in.preds[-i]) == 0, "no fixed predictors", paste(in.preds[-i], collapse = " + ")) cat(sep = "\n", paste("Run ", run, " (above) with ", above.preds, " is worse than Run ", best.run, ifelse(above.preds == "no fixed predictors", "", " also"), " including ", in.preds[i], ", ", p.flag, signif(compare[i], decimals), sep = ""))} run <- run + 1} repeat{ if (all(compare %in% c(-888, -999))){ if(any(in.preds %in% interactions)) best.model <- mod[[xxx <- (which(in.preds %in% interactions)[1])]] else best.model <- mod[[xxx <- 1]] best.run <- mod.run[xxx] if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run last step or this step, dropping first free term ", in.preds[xxx], "...", sep = ""))} in.preds <- in.preds[-xxx] error.count <- error.count + 1 break} if (all(compare %in% c(888, 999))){ if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run this step, best model from last step is Run ", best.run, sep = ""))} keep.going <- F error.count <- error.count + 1 break} if (max(compare) < 0){ best.model <- mod[[yyy <- which(compare == max(compare))[1]]] best.run <- mod.run[yyy] if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run last step, best model this step is Run ", best.run, ", dropping ", in.preds[yyy], "...", sep = ""))} in.preds <- in.preds[-yyy] error.count <- error.count + 1 break} if ((zzz <- max(compare[compare < 888])) >= threshold){ best.model <- mod[[yyy <- which(compare == zzz)[1]]] best.run <- mod.run[yyy] if (!silent) { cat("\n") cat(sep = "\n", paste("Dropping ", in.preds[yyy], "...", sep = ""))} in.preds <- in.preds[-yyy] break} if (zzz < threshold){ if (!silent) { cat("\n") cat(sep = "\n", paste("All remaining predictors are significant, best model from last step is Run ", best.run, sep = ""))} keep.going <- F break} error <- c(success, ful)}} if (length(in.preds) == 0){ in.preds <- "no fixed predictors" compare <- 1234} compare2 <- signif(compare[order(compare)], decimals) compare2[compare2 < .Machine$double.xmin] <- "~0" down0 <- paste( ifelse(length(random.i) > 0, paste(c(random.i, random.s), " [random, not tested] and ", sep = "", collapse = ""), ""), paste(in.preds[order(compare)], " (", ifelse(compare == 1234, "", compare2), ")", sep = "", collapse = " + "), "\n[p-values dropping from full model]", sep = "") down1 <- gsub("\\(999\\)", "[error]", down0) down2 <- gsub("\\(-999\\)", "[error]", down1) down3 <- gsub("\\(888\\)", "[main effect, not tested]", down2) down4 <- gsub("\\(-888\\)", "[main effect, not tested]", down3) down <<- down4 down.preds <<- in.preds fixed <<- fix continuous <<- cont interactions <<- int model <<- list(best.model, down, contr) down.model <<- best.model if (!silent) { cat("\n") cat(sep = "\n", paste("BEST STEP-DOWN MODEL OF RESPONSE", response, "IS WITH PREDICTOR(S):", down)) if (verbose) {cat("\n"); print(extract(down.model), row.names = F)} if (error.count > 0) cat(sep = "\n", paste("THERE WERE ERRORS ON", error.count, "OF", steppe + 1, "STEPS!"))} down.model} step.one <- function(){ if (is.null(variables)) return(cat(sep = "\n", "Please choose variables first.")) load.variables() fix <- fixed cont <- continuous int <- interactions one.preds <<- c(fixed, continuous, interactions) in.preds <- c(fixed, continuous, interactions) one.model <<- fit.model() compare <- vector() mod <- list() if (length(in.preds) > 0) for (i in 1:length(in.preds)){ if (in.preds[i] %in% sapply(in.preds[in.preds %in% int], split.term)){ mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error" compare[i] <- get.p(one.model, mod[[i]]) next} selection <- in.preds[-i] continuous <<- cont[cont %in% selection] fixed <<- fix[fix %in% selection] interactions <<- int[int %in% selection] mod[[i]] <- fit.model() compare[i] <- get.p(one.model, mod[[i]])} compare2 <- signif(compare[order(compare)], decimals) compare2[compare2 < .Machine$double.xmin] <- "~0" one0 <- paste(ifelse(length(random.i) > 0, paste(c(random.i, random.s), " [random, not tested] and ", sep = "", collapse = ""), ""), paste(in.preds[order(compare)], " (", compare2, ")", sep = "", collapse = " + "), sep = "") one1 <- gsub("\\(999\\)", "[error]", one0) one2 <- gsub("\\(-999\\)", "[error]", one1) one3 <- gsub("\\(888\\)", "[main effect, not tested]", one2) one4 <- gsub("\\(-888\\)", "[main effect, not tested]", one3) one <<- one4 fixed <<- fix continuous <<- cont interactions <<- int if (!silent){ cat("\n") cat(sep = "\n", paste("ONE-LEVEL ANALYSIS OF RESPONSE", response, "WITH PREDICTOR(S):", one)) if (verbose) {cat("\n"); print(extract(one.model), row.names = F)}} model <<- list(one.model, one, contr) one.model} step.up <- function(){ if (is.null(variables)) return(cat(sep = "\n", "Please choose variables first.")) load.variables() fix <- fixed cont <- continuous int <- interactions out.preds <- c(fixed, continuous, interactions) in.preds <- "no fixed predictors" p.add <- vector() keep.going <- T steppe <- 0 run <- 1 fixed <<- 1 continuous <<- 1 interactions <<- 1 best.model <- fit.model() if ("try-error" %in% class(best.model)) error.count <- 1 else error.count <- 0 best.run <- 0 if (!silent){ cat("\n") cat(sep = "\n", "STEPPING UP...") cat("\n") cat(sep = "\n", "STEP 0 - Run 0 - model with no fixed predictors") cat("\n") if (verbose) print(extract(best.model), row.names = F)} while (keep.going & length(out.preds) > 0){ steppe <- steppe + 1 if (thresh == "A") threshold <<- 0.05/length(out.preds) else threshold <<- thresh compare <- vector() mod <- list() mod.run <- vector() if (!silent){ cat("\n") cat(sep = "\n", paste("STEP", steppe, "- adding to Run", best.run, "model with", paste(in.preds, collapse = " + ")))} for (i in 1:length(out.preds)){ if (out.preds[i] %in% int & !all(split.term(out.preds[i]) %in% in.preds)){ mod[[i]] <- "interaction"; class(mod[[i]]) <- "try-error" compare[i] <- get.p(best.model, mod[[i]]) next} cat("\n") cat(sep = "\n", paste("Trying with ", out.preds[i], "...", sep = "")) selection <- c(in.preds, out.preds[i]) continuous <<- cont[cont %in% selection] fixed <<- fix[fix %in% selection] interactions <<- int[int %in% selection] mod[[i]] <- fit.model() mod.run[i] <- run compare[i] <- get.p(best.model, mod[[i]]) if (!silent) {cat("\n") if (verbose) print(extract(mod[[i]]), row.names = F) above.preds <- ifelse(all(in.preds == "no fixed predictors"), out.preds[i], paste(paste(in.preds, collapse = " + "), out.preds[i], sep = " + ")) cat(sep = "\n", paste("Run ", run, " (above) with ", above.preds, " is better than Run ", best.run, " without ", out.preds[i], ", ", p.flag, signif(compare[i], decimals), sep = ""))} run <- run + 1} repeat{ if (all(compare %in% c(-888, -999))){ if (any(!out.preds %in% int)) best.model <- mod[[xxx <- (which(!out.preds %in% int)[1])]] else best.model <- mod[[xxx <- 1]] best.run <- mod.run[xxx] if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run last step or this step, adding first free term ", out.preds[xxx], "...", sep = ""))} in.preds <- in.preds[in.preds!= "no fixed predictors"] in.preds <- c(in.preds, out.preds[xxx]) out.preds <- out.preds[!out.preds %in% in.preds] error.count <- error.count + 1 break} if (all(compare %in% c(888, 999))){ if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run this step, best model from last step is Run ", best.run, sep = ""))} keep.going <- F error.count <- error.count + 1 break} if (max(compare) < 0){ best.model <- mod[[yyy <- which(compare == max(compare))[1]]] best.run <- mod.run[yyy] if (!silent) { cat("\n") cat(sep = "\n", paste("No successful run last step, best model this step is Run ", best.run, ", adding ", out.preds[yyy], "...", sep = ""))} in.preds <- c(in.preds, out.preds[yyy]) out.preds <- out.preds[!out.preds %in% in.preds] error.count <- error.count + 1 break} if (min(compare) < threshold){ best.model <- mod[[zzz <- (which(compare == min(compare))[1])]] best.run <- mod.run[zzz] if (!silent) { cat("\n") cat(sep = "\n", paste("Adding ", out.preds[zzz], "...", sep = ""))} in.preds <- in.preds[in.preds!= "no fixed predictors"] in.preds <- c(in.preds, out.preds[zzz]) out.preds <- out.preds[!out.preds %in% in.preds] p.add <- c(p.add, signif(min(compare), decimals)) break} if (min(compare) >= threshold){ if (!silent) { cat("\n") cat(sep = "\n", paste("No significant improvement this step, best model from last step is Run ", best.run, sep = ""))} keep.going <- F break} error <- c(success, ful)}} if (length(p.add) == 0) p.add <- 1234 p.add2 <- signif(p.add, decimals) p.add2[p.add2 < .Machine$double.xmin] <- "~0" up0 <- paste(ifelse(length(random.i) > 0, paste(c(random.i, random.s), " [random, not tested] and ", sep = "", collapse = ""), ""), paste(in.preds, " (", ifelse(p.add == 1234, "", p.add2), ")", sep = "", collapse = " + "), "\n[p-values building from null model]", sep = "") up1 <- gsub("\\(999\\)", "[error]", up0) up2 <- gsub("\\(-999\\)", "[error]", up1) up3 <- gsub("\\(888\\)", "[main effect, not tested]", up2) up4 <- gsub("\\(-888\\)", "[main effect, not tested]", up3) up <<- up4 fixed <<- fix continuous <<- cont interactions <<- int up.preds <<- in.preds model <<- list(best.model, up, contr) up.model <<- best.model if (!silent) {cat("\n") cat(sep = "\n", paste("BEST STEP-UP MODEL OF RESPONSE", response, "IS WITH PREDICTOR(S):", up)) if (verbose) {cat("\n"); print(extract(up.model), row.names = F)} if (error.count > 0) cat(sep = "\n", paste("THERE WERE ERRORS ON", error.count, "OF", steppe+1, "STEPS!"))} up.model} step.updown <- function(){ if (is.null(variables)) return(cat(sep = "\n", "Please choose variables first.")) step.up() step.down() cat("\n") if (!silent){ cat(sep = "\n", paste("BEST STEP-UP MODEL OF RESPONSE", response, "IS WITH PREDICTOR(S):", up)) cat("\n") cat(sep = "\n", paste("BEST STEP-DOWN MODEL OF RESPONSE", response, "IS WITH PREDICTOR(S):", down)) cat("\n") if (all(up.preds %in% down.preds) & all(down.preds %in% up.preds)) cat(sep = "\n", "STEP-UP AND STEP-DOWN MATCH!\n") else cat(sep = "\n", "STEP-UP AND STEP-DOWN MISMATCH!!\n")}} sumify <- function(data, ncol){ data <- matrix(data, ncol = ncol) data[1, max(which(is.na(data[1, ])))] <- -sum(data[1, ], na.rm = T) if (nrow(data) > 2){ for (r in 2:(nrow(data) - 1)){ data[r, max(which(is.na(data[r, ])))] <- -sum(data[r, ], na.rm = T)}} if (nrow(data) > 1){ data[nrow(data), ] <- -colSums(data, na.rm = T)} as.vector(data)} test.chisq <- function(){ cat(sep = "\n", paste("Enter first deviance or log likelihood.")) repeat { d1 <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (!"try-error" %in% class(d1)) break cat(sep = "\n", paste("Invalid entry, please try again."))} cat(sep = "\n", paste("Enter second deviance or log likelihood.")) repeat { d2 <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (!"try-error" %in% class(d2)) break cat(sep = "\n", paste("Invalid entry, please try again."))} cat(sep = "\n", paste("If these were log likelihood values, press 1. Press Enter if they were deviances.")) repeat { ll.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(ll.pick) == 0) break if (ll.pick == 1) {d1 <-2*d1; d2 <-2*d2; break} cat(sep = "\n", paste("Invalid entry, please try again."))} cat(sep = "\n", paste("Enter difference in degrees of freedom.")) repeat { df <- try(as.integer(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (!"try-error" %in% class(df)) break cat(sep = "\n", paste("Invalid entry, please try again."))} xx <- abs(d1-d2) cat(sep = "\n", paste("Chi-square = ", xx, ", df = ", df, ", p = ", signif(pchisq(xx, df, lower = F), decimals), sep = ""))} transform <- function(){ cat(sep = "\n", paste("Continuous variable(s) to center or transform? (press Enter to exit)", get.menu(names(datafile)), sep = "")) repeat {transform.pick <- try(as.numeric(scan(what = "character", nmax = length(names(datafile)), quiet = T)), silent = T) transformees <- names(datafile)[transform.pick] if (!NA %in% transformees) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(transformees) > 0) for (i in transformees){ if (!class(datafile[, i]) %in% c("integer", "numeric")){ cat(sep = "\n", paste("Please change class of", i, "to 'continuous' before transforming.")) next} cat(sep = "\n", paste("Do what to ", i, "? (1-log 2-logit 3-center (rezero) 4-z-score Enter-nothing)", sep = "")) repeat {trans.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(trans.pick) == 0) break if (trans.pick %in% 1:6) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(trans.pick) == 0) break if (trans.pick == 1){ trans.col <- paste("log", i, sep = ".") if (any(datafile[, i] <= 0)){ cat(sep = "\n", "Not so fast! This variable has values that are zero or negative.") break} datafile[, trans.col] <<- log(datafile[, i])} if (trans.pick == 2){ trans.col <- paste("logit", i, sep = ".") if (max(datafile[, i], na.rm = T) <= 1 & min(datafile[, i], na.rm = T) >= 0){ datafile[, trans.col] <<- qlogis(datafile[, i])} else datafile[, trans.col] <<- NA} if (trans.pick == 3){ trans.col <- paste("ctr", i, sep = ".") cat(sep = "\n", paste("Center (rezero) on what value? (Enter - mean (", signif(center.mean <- mean(datafile[, i], na.rm = T), 5), "))", sep = "")) repeat {center.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (!is.na(center.mean) & length(center.pick) == 0) {center.pick <- center.mean; break} if (is.numeric(center.pick) & length(center.pick) > 0) break cat(sep = "\n", "Invalid entry, please try again.")} datafile[, trans.col] <<- datafile[, i]-center.pick} if (trans.pick == 4){ trans.col <- paste("z", i, sep = ".") datafile[, trans.col] <<- (datafile[, i]-mean(datafile[, i], na.rm = T))/ sd(datafile[, i], na.rm = T)}}} trim <- function(){ if (is.null(model)) return(cat(sep = "\n", "No model loaded.")) orig.datafile <- datafile orig.model <- model ff <- vector() if (length(fixed) == 0) fixed3 <- character() if (length(fixed) > 0) for (f in fixed){ ff <- c(ff, length(levels(datafile[, f]))) fixed3 <- fixed[which(ff > 2)]} if (length(fixed3) == 0){ cat("\n") cat(sep = "\n", "This model has no fixed-effect predictors (with 3 or more levels) to trim.") return()} big.model <- step.one() if (length(fixed3) == 1){ fixeds <- fixed3} if (length(fixed3) > 1){ cat("\n") cat(sep = "\n", paste("Choose predictors to combine levels within; press Enter for all", get.menu(fixed3), ".", sep = "")) repeat {fixed.pick <- try(as.numeric(scan(what = "character", nmax = length(fixed3), quiet = T)), silent = T) if (length(fixed.pick) == 0) {fixeds <- fixed3; break} if (!NA %in% fixed3[fixed.pick]) {fixeds <- fixed3[fixed.pick]; break} cat(sep = "\n", "Invalid entry, please try again.")}} chooselevels <- vector() for (f in fixeds){ for (l in levels(datafile[, f])){ chooselevels <- c(chooselevels, paste(f, l, sep = "."))}} cat("\n") cat(sep = "\n", paste("Choose levels to consider combining; press Enter for all", get.menu(chooselevels), ".", sep = "")) repeat {levels.pick <- try(as.numeric(scan(what = "character", nmax = length(chooselevels), quiet = T)), silent = T) if (length(levels.pick) == 0) break if (length(levels.pick) == 1) {cat(sep = "\n", "Choose at least two levels."); next} if (!NA %in% chooselevels[levels.pick]) {chooselevels <- chooselevels[levels.pick]; break} cat(sep = "\n", "Invalid entry, please try again.")} for (loop in 1:1000){ ff <- vector() for (f in fixeds) ff <- c(ff, length(levels(datafile[, f]))) fixeds <- fixeds[which(ff > 2)] if (length(fixeds) == 0){ cat("\n") cat(sep = "\n", "This model has no fixed-effect predictors (with 3 or more levels) to trim.") break} pairframe <- data.frame() choices <- data.frame() for (f in fixeds){ l <- levels(datafile[, f]) choices <- rbind(choices, cbind(left = combn(l, 2)[1, ], right = combn(l, 2)[2, ])) for (ch in 1:length(choices[, 1])){ if (paste(f, choices$left[ch], sep = ".") %in% chooselevels & paste(f, choices$right[ch], sep = ".") %in% chooselevels){ pairframe <- rbind(pairframe, cbind(factor = f, left = as.character(choices$left[ch]), right = as.character(choices$right[ch])))}}} if (length(pairframe) == 0) break for (i in 1:length(pairframe[, 1])){ little.datafile <- datafile levs <- levels(little.datafile[, as.character(pairframe$factor[i])]) leftnum <- which(levs == pairframe$left[i]) rightnum <- which(levs == pairframe$right[i]) levs[c(leftnum, rightnum)] <- paste(pairframe$left[i], pairframe$right[i], sep = "+") levels(little.datafile[, as.character(pairframe$factor[i])]) <- levs pairframe$p[i] <- signif(get.p(big.model, fit.model(little.datafile)), decimals)} pairframe <- cbind(row = 1:length(pairframe[, 1]), pairframe) pairframe2 <- cbind(pairframe, " ", pairframe[, -1], pairframe[, 1]) names(pairframe2) <- c(names(pairframe), " ", paste(" ", c(names(pairframe)[-1], names(pairframe[1])), sep = "")) names(pairframe2)[c(2, 7)] <- "predictor" for (pf in c(7, 8, 9, 11, 10)) pairframe2[, pf] <- pairframe2[order(pairframe2[, 10]), pf] cat("\n") print(pairframe2) cat("\n") cat(sep = "\n", "Type a row number to combine a pair of levels ('left' and 'right'). Press Enter to stop.") repeat {combine.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(combine.pick) == 0) break if (combine.pick %in% pairframe$row) break cat(sep = "\n", "Invalid entry, please try again.")} if (length(combine.pick) == 0) break pair <- pairframe[combine.pick, ] levs <- levels(datafile[, as.character(pair$factor)]) leftnum <- which(levs == pair$left) rightnum <- which(levs == pair$right) cat("\n") cat(sep = "\n", paste("Type new name for combination or press Enter to use ", (lr <- paste(pair$left, pair$right, sep = "+")), "?", sep = "")) repeat {combo.pick <- try(scan(what = "character", nmax = 1, quiet = T), silent = T) if (length(combo.pick) == 0){ levs[c(leftnum, rightnum)] <- lr chooselevels <- c(chooselevels, paste(pair$factor, lr, sep = ".")) break} else {tlr <- try(levs[c(leftnum, rightnum)] <- combo.pick) if (!"try-error" %in% class(tlr)){ chooselevels <- c(chooselevels, paste(pair$factor, combo.pick, sep = ".")) break}} cat(sep = "\n", "Invalid entry, please try again.")} levels(datafile[, as.character(pair$factor)]) <<- levs big.model <- step.one()} cat("\n") cat(sep = "\n", "Press 1 to keep trimmed data file and model, press Enter to keep originals.") repeat {keep.pick <- try(as.numeric(scan(what = "character", nmax = 1, quiet = T)), silent = T) if (length(keep.pick) == 0){ datafile <<- orig.datafile model <<- orig.model break} if (keep.pick == 1){ break} cat(sep = "\n", "Invalid entry, please try again.")}} vif.dej <- function(mod){ ### hacked from car::vif to handle mixed models ### term-wise results, not parameter-wise as in vif.mer if ("glm" %in% class(mod)) coefs <- coef if (class(mod)[1] %in% c("lmerMod", "glmerMod")) coefs <- fixef if (any(is.na(coefs(mod)))) stop("error: NA's") v <- vcov(mod) assign <- attr(model.matrix(mod), "assign") if (names(coefs(mod)[1]) == "(Intercept)"){ v <- v[-1, -1] assign <- assign[-1]} else stop("no intercept") terms <- labels(terms(mod)) n.terms <- length(terms) if (n.terms < 2) stop("model contains fewer than 2 terms") R <- cov2cor(v) detR <- det(R) result <- vector(length = n.terms) names(result) <- terms for (term in 1:n.terms){ subs <- which(assign == term) result[term] <- det(as.matrix(R[subs, subs])) * det(as.matrix(R[-subs, -subs])) / detR} result} ### following is copied from MuMIn package because MuMIn package is not installing properly r.squaredGLMM <- function(x){ if ("glm" %in% class(x)){ fam <- family(x) return(rsqGLMM(fam, varFx = var(as.vector(model.matrix(x) %*% coef(x))), varRe = 0, varResid = sum(if (is.null(x$weights)) resid(x)^2 else resid(x)^2 * x$weights)/df.residual(x), beta0 = if (fam$family == "poisson" && fam$link == "log") log(mean(model.response(model.frame(x)))) else NULL)) } if (class(x)[1] %in% c("lmerMod", "glmerMod")){ fam <- family(x) useObsLevVar <- (fam$family == "poisson" && fam$link == "log") || fam$family == "binomial" if (useObsLevVar && !any(sapply(x@flist, nlevels) == nobs(x))) { cl <- get_call(x) frm <- formula(x) nRowData <- eval(call("eval", as.expression(call("NROW", cl$formula[[2L]])), envir = cl$data), envir = environment(frm), enclos = parent.frame()) fl <- length(frm) frx <- . ~ . + 1 frx[[3L]][[3L]] <- call("(", call("|", 1, call("gl", nRowData, 1))) cl$formula <- update.formula(frm, frx) x <- tryCatch(eval(cl, envir = environment(frm), enclos = parent.frame()), error = function(e) { cry(conditionCall(e), conditionMessage(e), warn = TRUE) cry(cl, "fitting model with the observation-level random effect term failed. Add the term manually") }) message("The result is correct only if all data used by the model ", "has not changed since model was fitted.", domain = "R-MuMIn") } mmAll <- model.matrix(ranform(formula(x)), data = model.frame(x)) vc <- lme4::VarCorr(x) n <- nrow(mmAll) fx <- lme4::fixef(x) fxpred <- as.vector(model.matrix(x) %*% fx) if (useObsLevVar) { vname <- names(x@flist)[sapply(x@flist, nlevels) == n][1L] varResid <- vc[[vname]][1L] beta0 <- mean(fxpred) vc <- vc[names(vc) != vname] } else { varResid <- attr(vc, "sc")^2 beta0 <- NULL } if (!all(unlist(sapply(vc, rownames), use.names = FALSE) %in% colnames(mmAll))) stop("random term names do not match those in model matrix. \n", "Have 'options(contrasts)' changed since the model was fitted?") varRe <- ifelse(length(vc) == 0L, 0L, sum(sapply(vc, function(sig){ mm1 <- mmAll[, rownames(sig), drop = FALSE] sum(matmultdiag(mm1 %*% sig, ty = mm1))/n}))) return(rsqGLMM(fam = family(x), varFx = var(fxpred), varRe = varRe, varResid = varResid, beta0 = beta0)) }} rsqGLMM <- function (fam, varFx, varRe, varResid, beta0){ varDistr <- switch(paste(fam$family, fam$link, sep = "."), gaussian.identity = 0, binomial.logit = 3.28986813369645, binomial.probit = 1, poisson.log = { expBeta0 <- exp(beta0) if (expBeta0 < 6) cry(sys.call(-1L), "exp(beta0) of %0.1f is too close to zero, estimate may be unreliable \n", expBeta0, warn = TRUE) log1p(1/expBeta0) }, poisson.sqrt = 0.25, cry(sys.call(-1L), "do not know how to calculate variance for this family/link combination")) varTot <- sum(varFx, varRe) res <- c(varFx, varTot)/(varTot + varDistr + varResid) names(res) <- c("R2m", "R2c") res } ranform <- function(form){ ans <- update.formula(reformulate(vapply(lapply(findbars(form), "[[", 2L), deparse, "")), ~. + 1) environment(ans) <- environment(form) ans } findbars <- function (term){ if (is.name(term) || !is.language(term)) return(NULL) if (term[[1L]] == as.name("(")) return(findbars(term[[2L]])) if (!is.call(term)) stop("term must be of class call") if (term[[1L]] == as.name("|")) return(term) if (length(term) == 2L) return(findbars(term[[2L]])) c(findbars(term[[2L]]), findbars(term[[3L]])) } matmultdiag <- function (x, y, ty = t(y)){ if (ncol(x) != ncol(ty)) stop("non-conformable arguments") if (nrow(x) != nrow(ty)) stop("result is not a square matrix") return(rowSums(x * ty)) } AICc <- function (object, ..., k = 2, REML = NULL){ aic(list(object, ...), 1, k, REML, function(ll, chat, k, npar.adj) { no <- get.nobs(ll, error = FALSE) df <- attr(ll, "df") ret <- (-2 * c(ll)) + (k * df) * (1 + ((df + 1)/(no - df - 1))) return(ret) }, "AICc") } aic <- function (objectlist, chat, k, REML, ICFun, ICName){ if (chat < 1) { warning("'chat' given is < 1, increased to 1") chat <- 1 } npar.adj <- if (chat == 1) 0 else 1 llCall <- call("logLik", as.name("object")) if (!is.null(REML)) llCall$REML <- REML ll <- function(object) fixLogLik(NA, object) body(ll)[[2L]] <- llCall if (length(objectlist) > 1L) { lls <- lapply(objectlist, ll) val <- data.frame(df = vapply(lls, attr, 1, "df"), ic = sapply(lls, ICFun, chat = chat, k = k, npar.adj = npar.adj)) Call <- match.call(sys.function(sys.parent()), call = sys.call(sys.parent())) Call$chat <- Call$REML <- Call$k <- NULL dimnames(val) <- list(as.character(Call[-1L]), c("df", ICName)) return(val) } else { return(ICFun(ll(objectlist[[1L]]), chat = chat, k = k, npar.adj = npar.adj)) } } get.nobs <- function (llik, error = FALSE){ no <- attr(llik, "nall") if (is.null(no)) no <- attr(llik, "nobs") if (error && is.null(no)) stop("'logLik' object must have a \"nobs\" attribute") no } fixLogLik <- function (ll, object){ if (is.null(attr(ll, "nall")) && is.null(attr(ll, "nobs"))) attr(ll, "nobs") <- nobs(object) ll } get_call <- function (x) { rval <- if (isS4(x)) { if (any(i <- (sln <- c("call", "CALL", "Call")) %in% methods::slotNames(x))) slot(x, sln[i][1L]) else if (!is.null(attr(x, "call"))) attr(x, "call") else NULL } else { if (!is.atomic(x) && (i <- match("call", names(x), nomatch = 0L)) != 0L) { x[[i]] } else if (!is.null(attr(x, "call"))) { attr(x, "call") } else NULL } if (is.null(rval)) stats::getCall(x) else rval } ##### end of section copied from MuMIn package #####