diff --git a/R/commonglm.R b/R/commonglm.R index abc1f11d..fb4a7321 100644 --- a/R/commonglm.R +++ b/R/commonglm.R @@ -5,7 +5,7 @@ type <- "binomial" if (type == "binomial") { # Logistic regression - ff <- .createGlmFormula(options) + ff <- .createLogisticRegFormula(options) nf <- .createNullFormula(options) # calculate null and full models nullMod <- glm(nf, family = "binomial", data = dataset) @@ -18,7 +18,53 @@ return(glmRes) } -.createGlmFormula <- function(options) { +# create GLM formula +.createGLMFormula <- function(options, nullModel = FALSE) { + + leftTerm <- options$dependent + modelTerms <- options$modelTerms + includeIntercept <- options$interceptTerm + offsetTerm <- options$offset + + if (includeIntercept) + rightTerms <- "1" + else + rightTerms <- "-1" + + if (offsetTerm != "") + rightTerms <- c(rightTerms, paste("offset(", offsetTerm, ")", sep = "")) + + if (length(modelTerms) == 0) { + f <- formula(paste(leftTerm, "~", rightTerms)) + } else { + + if (nullModel) { + for (i in seq_along(modelTerms)) { + nuisance <- modelTerms[[i]][["isNuisance"]] + if (!is.null(nuisance) && nuisance) { + term <- modelTerms[[i]][["components"]] + if (length(term) == 1) + rightTerms <- c(rightTerms, term) + else + rightTerms <- c(rightTerms, paste(term, collapse = ":")) + } + } + + } else { + for (i in seq_along(modelTerms)) { + term <- modelTerms[[i]][["components"]] + if (length(term) == 1) + rightTerms <- c(rightTerms, term) + else + rightTerms <- c(rightTerms, paste(term, collapse = ":")) + } + } + f <- formula(paste(leftTerm, "~", paste(rightTerms, collapse = "+"))) + } + return(f) +} + +.createLogisticRegFormula <- function(options) { # this function outputs a formula name with base64 values as varnames f <- NULL @@ -123,6 +169,312 @@ return(modlist) } +# as explained in ?is.integer +.is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol + + +.glmComputeModel <- function(jaspResults, dataset, options) { + if (!is.null(jaspResults[["glmModels"]])) + return(jaspResults[["glmModels"]]$object) + + # make formulas for the full model and the null model + ff <- .createGLMFormula(options, nullModel = FALSE) + environment(ff) <- environment() + nf <- .createGLMFormula(options, nullModel = TRUE) + environment(nf) <- environment() + + weights <- dataset[[options[["weights"]]]] + + if (options$family != "other") { + # specify family and link + if (options$family == "bernoulli") { + family <- "binomial" + } + else if (options$family == "gamma") { + family <- "Gamma" + } + else if (options$family == "inverseGaussian") { + family <- "inverse.gaussian" + } + else { + family <- options$family + } + + familyLink <- eval(call(family, link = options$link)) + # compute full and null models + fullModel <- try(stats::glm(ff, + family = familyLink, + data = dataset, + weights = weights)) + nullModel <- try(stats::glm(nf, + family = familyLink, + data = dataset, + weights = weights)) + } else { + + if (options$otherGlmModel == "multinomialLogistic") { + fullModel <- try(VGAM::vglm(ff, + family = VGAM::multinomial(), + data = dataset, + weights = weights)) + nullModel <- try(VGAM::vglm(nf, + family = VGAM::multinomial(), + data = dataset, + weights = weights)) + } + + if (options$otherGlmModel == "ordinalLogistic") { + fullModel <- try(VGAM::vglm(ff, + family = VGAM::cumulative(link = "logitlink", parallel = TRUE), + data = dataset, + weights = weights)) + nullModel <- try(VGAM::vglm(nf, + family = VGAM::cumulative(link = "logitlink", parallel = TRUE), + data = dataset, + weights = weights)) + } + + if (options$otherGlmModel == "firthLogistic") { + fullModel <- try(logistf::logistf(ff, + data = dataset, + pl = TRUE, + firth = TRUE, + weights = weights)) + nullModel <- try(logistf::logistf(nf, + data = dataset, + pl = TRUE, + firth = TRUE, + weights = weights)) + } + } + + if (jaspBase::isTryError(fullModel)) { + msg <- .glmFitErrorMessageHelper(fullModel) + msg <- gettextf("The full model could not be fitted to the data, with the following error message: '%s'", msg) + jaspBase::.quitAnalysis(msg) + } + + if (jaspBase::isTryError(nullModel)) { + msg <- .glmFitErrorMessageHelper(nullModel) + msg <- gettextf("The null model could not be fitted to the data, with the following error message: '%s'", msg) + jaspBase::.quitAnalysis(msg) + } + + glmModels <- list("nullModel" = nullModel, + "fullModel" = fullModel) + # combine both models + jaspResults[["glmModels"]] <- createJaspState() + jaspResults[["glmModels"]]$dependOn(optionsFromObject = jaspResults[["modelSummary"]]) + jaspResults[["glmModels"]]$object <- glmModels + + return(glmModels) +} + +.glmFitErrorMessageHelper <- function(tryError) { + msg <- jaspBase::.extractErrorMessage(tryError) + if(msg == "cannot find valid starting values: please specify some") + msg <- gettext("Could not find valid starting values. Check feasibility of the model to fit the data.") + + return(msg) +} + + +.hasNuisance <- function(options) { + return(any(sapply(options$modelTerms, function(x) x[["isNuisance"]]))) +} + +.glmCreatePlotPlaceholder <- function(container, index, title) { + jaspPlot <- createJaspPlot(title = title) + jaspPlot$status <- "running" + container[[index]] <- jaspPlot + return() +} + + +.glmStdResidCompute <- function(model, residType, options) { + if (residType == "deviance") { + stdResid <- rstandard(model) + } + else if (residType == "Pearson") { + phiEst <- summary(model)[["dispersion"]] + resid <- resid(model, type="pearson") + stdResid <- resid / sqrt(phiEst * (1 - hatvalues(model))) + } + else if (residType == "quantile") { + jaspBase::.setSeedJASP(options) + resid <- statmod::qresid(model) + stdResid <- resid / sqrt(1 - hatvalues(model)) + } else { + stdResid <- NULL + } + + return(stdResid) +} + +.glmInsertPlot <- function(jaspPlot, func, ...) { + p <- try(func(...)) + + if (inherits(p, "try-error")) { + errorMessage <- .extractErrorMessage(p) + jaspPlot$setError(gettextf("Plotting is not possible: %s", errorMessage)) + } else { + jaspPlot$plotObject <- p + jaspPlot$status <- "complete" + } +} + +# Plots: Residuals Q-Q +.glmPlotResQQ <- function(jaspResults, dataset, options, ready, position) { + + plotNames <- c("devianceResidualQqPlot", "pearsonResidualQqPlot", "quantileResidualQqPlot") + if (!ready || !any(unlist(options[plotNames]))) + return() + + residNames <- c("deviance", "Pearson", "quantile") + + glmPlotResQQContainer <- createJaspContainer(gettext("Q-Q Plots")) + glmPlotResQQContainer$dependOn(optionsFromObject = jaspResults[["modelSummary"]], + options = c(plotNames, "seed", "setSeed")) + glmPlotResQQContainer$position <- position + jaspResults[["diagnosticsContainer"]][["glmPlotResQQ"]] <- glmPlotResQQContainer + + + if (!is.null(jaspResults[["glmModels"]])) { + glmFullModel <- jaspResults[["glmModels"]][["object"]][["fullModel"]] + for (i in 1:length(plotNames)) { + if (options[[plotNames[[i]]]]) { + .glmCreatePlotPlaceholder(glmPlotResQQContainer, + index = plotNames[[i]], + title = gettextf("Q-Q plot: Standardized %1s residuals", residNames[[i]])) + + .glmInsertPlot(glmPlotResQQContainer[[plotNames[[i]]]], + .glmFillPlotResQQ, + residType = residNames[[i]], + model = glmFullModel, + options = options) + } + } + } + return() +} + +.glmFillPlotResQQ <- function(residType, model, options) { + + # compute residuals + stdResid <- .glmStdResidCompute(model = model, residType = residType, options = options) + + p <- jaspGraphs::plotQQnorm(stdResid, ablineColor = "darkred", ablineOrigin = TRUE, identicalAxes = TRUE) + + return(p) +} + + +.regressionExportResiduals <- function(container, model, dataset, options, ready) { + + if (isFALSE(options[["residualsSavedToData"]])) + return() + + if (is.null(container[["residualsSavedToDataColumn"]]) && options[["residualsSavedToDataColumn"]] != "") { + + residuals <- model[["residuals"]] # extract residuals + + container[["residualsSavedToDataColumn"]] <- createJaspColumn(columnName = options[["residualsSavedToDataColumn"]]) + container[["residualsSavedToDataColumn"]]$dependOn(options = c("residualsSavedToDataColumn", "residualsSavedToData")) + container[["residualsSavedToDataColumn"]]$setScale(residuals) + + } + +} + +.constInfoTransform <- function(family, x) { + switch(family, + "bernoulli" = 1/(sin(sqrt(x))), + "binomial" = 1/(sin(sqrt(x))), + "poisson" = sqrt(x), + "Gamma" = log(x), + "inverse.gaussian" = 1/sqrt(x), + "gaussian" = x) +} + +.constInfoTransName <- function(family) { + switch(family, + "bernoulli" = expression(sin^-1 * sqrt("Fitted values")), + "binomial" = expression(sin^-1 * sqrt("Fitted values")), + "poisson" = expression(sqrt("Fitted values")), + "Gamma" = expression(log("Fitted values")), + "inverse.gaussian" = expression(1/sqrt("Fitted Values")), + "gaussian" = "Fitted values") +} + +.capitalize <- function(x) { + x_new <- paste(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)), sep="") + return(x_new) +} + +# function for multicollineary statistics, taken from the source code of car::vif +# car version: 3.0-13 see https://rdrr.io/cran/car/src/R/vif.R +# the reason is to remove the dependency on the car package. When importing the entire car package, +# JASP crashes when loading the regression module (at least on Windows 10). + +.vif.default <- function(mod, ...) { + # modified to fix bug: https://github.com/jasp-stats/jasp-test-release/issues/2487 + if (length(coef(mod)) == 0) { + stop(gettext("There are no predictors to test for multicollinearity")) + } + # end modification + if (any(is.na(coef(mod)))) + stop ("there are aliased coefficients in the model") + v <- vcov(mod) + assign <- attr(model.matrix(mod), "assign") + if (names(coefficients(mod)[1]) == "(Intercept)") { + v <- v[-1, -1] + assign <- assign[-1] + } + else warning("No intercept: vifs may not be sensible.") + 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 <- matrix(0, n.terms, 3) + rownames(result) <- terms + colnames(result) <- c("GVIF", "Df", "GVIF^(1/(2*Df))") + for (term in 1:n.terms) { + subs <- which(assign == term) + result[term, 1] <- det(as.matrix(R[subs, subs])) * + det(as.matrix(R[-subs, -subs])) / detR + result[term, 2] <- length(subs) + } + if (all(result[, 2] == 1)) result <- result[, 1] + else result[, 3] <- result[, 1]^(1/(2 * result[, 2])) + result +} + +# message for estimated marginal means table +.emmMessageTestNull <- function(value) gettextf("P-values correspond to test of null hypothesis against %s.", value) +.emmMessageAveragedOver <- function(terms) gettextf("Results are averaged over the levels of: %s.",paste(terms, collapse = ", ")) +.messagePvalAdjustment <- function(adjustment) { + if (adjustment == "none") { + return(gettext("P-values are not adjusted.")) + } + adjustment <- switch(adjustment, + "holm" = gettext("Holm"), + "hommel" = gettext("Homel"), + "hochberg" = gettext("Hochberg"), + "mvt" = gettext("Multivariate-t"), + "tukey" = gettext("Tukey"), + "BH" = gettext("Benjamini-Hochberg"), + "BY" = gettext("Benjamini-Yekutieli"), + "scheffe" = gettext("Scheffé"), + "sidak" = gettext("Sidak"), + "dunnettx" = gettext("Dunnett"), + "bonferroni" = gettext("Bonferroni") + ) + return(gettextf("P-values are adjusted using %s adjustment.", adjustment)) +} + + + .confusionMatAddColInfo <- function(table, levs, type) { table$addColumnInfo(name = "pred0", title = paste0(levs[1]), type = type, overtitle = gettext("Predicted")) table$addColumnInfo(name = "pred1", title = paste0(levs[2]), type = type, overtitle = gettext("Predicted")) @@ -383,6 +735,158 @@ return(robustSE) } + +# Table: Influential cases +.glmInfluenceTable <- function(jaspResults, model, dataset, options, ready, position, linRegAnalysis = FALSE) { + + tableOptionsOn <- c(options[["dfbetas"]], + options[["dffits"]], + options[["covarianceRatio"]], + options[["leverage"]], + options[["mahalanobis"]]) + + nModels <- length(options$modelTerms) + if (!ready || !options[["residualCasewiseDiagnostic"]] || + length(unlist(options$modelTerms[[nModels]][["components"]])) == 0) + return() + + + tableOptions <- c("dfbetas", "dffits", "covarianceRatio", "leverage", "mahalanobis") + tableOptionsClicked <- tableOptions[tableOptionsOn] + tableOptionsClicked <- c("cooksDistance", tableOptionsClicked) + + if (is.null(jaspResults[["influenceTable"]])) { + influenceTable <- createJaspTable(gettext("Table: Influential Cases")) + influenceTable$dependOn(optionsFromObject = jaspResults[["modelSummary"]], + options = tableOptions) + influenceTable$dependOn(c("residualCasewiseDiagnostic", "residualCasewiseDiagnosticType", + "residualCasewiseDiagnosticZThreshold", "residualCasewiseDiagnosticCooksDistanceThreshold")) + influenceTable$position <- position + influenceTable$showSpecifiedColumnsOnly <- TRUE + jaspResults[["influenceTable"]] <- influenceTable + } + + tableOptionToColName <- function(x) { + switch(x, + "dfbetas" = "DFBETAS", + "dffits" = "DFFITS", + "covarianceRatio" = "Covariance Ratio", + "cooksDistance" = "Cook's Distance", + "leverage" = "Leverage", + "mahalanobis" = "Mahalanobis") + } + + if (is.null(model)) { + for (option in tableOptionsClicked) { + colTitle <- tableOptionToColName(option) + influenceTable$addColumnInfo(name = option, title = gettext(colTitle), type = "number") + } + } else { + + colNameList <- c() + influenceTable$addColumnInfo(name = "caseN", title = "Case Number", type = "integer") + influenceTable$addColumnInfo(name = "stdResidual", title = gettext("Std. Residual"), type = "number", format = "dp:3") + influenceTable$addColumnInfo(name = "dependent", title = options$dependent, type = "number") + influenceTable$addColumnInfo(name = "predicted", title = gettext("Predicted Value"), type = "number") + influenceTable$addColumnInfo(name = "residual", title = gettext("Residual"), type = "number", format = "dp:3") + + alwaysPresent <- c("caseN", "stdResidual", "dependent", "predicted", "residual") + for (option in tableOptionsClicked) { + if (option == "dfbetas") { + predictors <- names(model$coefficients) + for (predictor in predictors) { + dfbetasName <- gettextf("DFBETAS_%1s", predictor) + colNameList <- c(colNameList, dfbetasName) + if (predictor == "(Intercept)") + dfbetasTitle <- gettext("DFBETAS:Intercept") + else + dfbetasTitle <- gettextf("DFBETAS:%1s", gsub(":", "*", predictor)) + influenceTable$addColumnInfo(name = dfbetasName, title = dfbetasTitle, type = "number") + } + } else { + colNameList <- c(colNameList, option) + colTitle <- tableOptionToColName(option) + influenceTable$addColumnInfo(name = option, title = gettext(colTitle), type = "number") + } + } + .glmInfluenceTableFill(influenceTable, dataset, options, ready, + model = model, + influenceMeasures = tableOptionsClicked, + colNames = c(colNameList, alwaysPresent)) + } +} + +.glmInfluenceTableFill <- function(influenceTable, dataset, options, ready, model, influenceMeasures, colNames) { + + influenceRes <- influence.measures(model) + nDFBETAS <- length(names(model$coefficients)) + + optionToColInd <- function(x, nDFBETAS) { + switch(x, + "dfbetas" = 1:nDFBETAS, + "dffits" = (nDFBETAS+1), + "covarianceRatio" = (nDFBETAS+2), + "cooksDistance" = (nDFBETAS+3), + "leverage" = (nDFBETAS+4))} + + colInd <- c() + for (measure in influenceMeasures) { + colInd <- c(colInd, optionToColInd(measure, nDFBETAS)) + } + + tempResult <- influenceRes[["infmat"]][, colInd] + resultContainsNaN <- any(is.nan(tempResult)) + tempResult[which(is.nan(tempResult))] <- NA + influenceResData <- as.data.frame(tempResult) + colnames(influenceResData)[1:length(colInd)] <- colNames[1:length(colInd)] + + influenceResData[["caseN"]] <- seq.int(nrow(influenceResData)) + influenceResData[["stdResidual"]] <- rstandard(model) + influenceResData[["dependent"]] <- model.frame(model)[[options$dependent]] + influenceResData[["predicted"]] <- model$fitted.values + influenceResData[["residual"]] <- model$residual + # browser() + modelMatrix <- as.data.frame(model.matrix(model)) + modelMatrix <- modelMatrix[colnames(modelMatrix) != "(Intercept)"] + influenceResData[["mahalanobis"]] <- mahalanobis(modelMatrix, center = colMeans(modelMatrix), cov = cov(modelMatrix)) + + if (options$residualCasewiseDiagnosticType == "cooksDistance") + index <- which(abs(influenceResData[["cooksDistance"]]) > options$residualCasewiseDiagnosticCooksDistanceThreshold) + else if (options$residualCasewiseDiagnosticType == "outliersOutside") + index <- which(abs(influenceResData[["stdResidual"]]) > options$residualCasewiseDiagnosticZThreshold) + else # all + index <- seq.int(nrow(influenceResData)) + + # funky statement to ensure a df even if only 1 row + influenceResSig <- subset(influenceRes[["is.inf"]], 1:nrow(influenceResData) %in% index, select = colInd) + colnames(influenceResSig) <- colNames[1:length(colInd)] + influenceResData <- influenceResData[index, ] + + if (length(index) == 0) + influenceTable$addFootnote(gettext("No influential cases found.")) + else { + influenceTable$setData(influenceResData) + # if any other metrix show influence, add footnotes: + if (sum(influenceResSig, na.rm = TRUE) > 0) { + for (thisCol in colnames(influenceResSig)) { + if (sum(influenceResSig[, thisCol], na.rm = TRUE) > 0) + influenceTable$addFootnote( + gettext("Potentially influential case, according to the selected influence measure."), + colNames = thisCol, + rowNames = rownames(influenceResData)[influenceResSig[, thisCol]], + symbol = "*" + ) + } + } + + if (resultContainsNaN) { + influenceTable$addFootnote( + gettext("Influence measures could not be computed for some cases due to extreme values, try another measure.") + ) + } + } +} + .casewiseDiagnosticsLogisticRegression <- function(dataset, model, options) { last <- length(model) diff --git a/R/correlation.R b/R/correlation.R index c2072099..df9bcdb3 100644 --- a/R/correlation.R +++ b/R/correlation.R @@ -115,7 +115,7 @@ CorrelationInternal <- function(jaspResults, dataset, options){ mainTable$dependOn(c("variables", "partialOutVariables", "pearson", "spearman", "kendallsTauB", "pairwiseDisplay", "significanceReport", "significanceFlagged", "sampleSize", - "ci", "ciLevel", + "ci", "ciLevel", "covariance", "vovkSellke", "alternative", "naAction", "ciBootstrap", "ciBootstrapSamples", "effectSize")) @@ -201,7 +201,14 @@ CorrelationInternal <- function(jaspResults, dataset, options){ if(options$significanceReport) mainTable$addColumnInfo(name = paste0(test, "_p.value"), title = gettext("p"), type = "pvalue", overtitle = overtitle) - + + if(options$vovkSellke){ + mainTable$addColumnInfo(name = paste0(test, "_vsmpr"), title = gettext("VS-MPR"), type = "number", overtitle = overtitle) + + mainTable$addFootnote(message = .corrGetTexts()$footnotes$VSMPR, symbol = "\u2020", colNames = paste0(test, "_vsmpr")) + mainTable$addCitation(.corrGetTexts()$references$Sellke_etal_2001) + } + if(options$ci){ mainTable$addColumnInfo(name = paste0(test, "_lower.ci"), title = gettextf("Lower %s%% CI", 100*options$ciLevel), type = "number", @@ -211,19 +218,16 @@ CorrelationInternal <- function(jaspResults, dataset, options){ overtitle = overtitle) } - if(options$vovkSellke){ - mainTable$addColumnInfo(name = paste0(test, "_vsmpr"), title = gettext("VS-MPR"), type = "number", overtitle = overtitle) - - mainTable$addFootnote(message = .corrGetTexts()$footnotes$VSMPR, symbol = "\u2020", colNames = paste0(test, "_vsmpr")) - mainTable$addCitation(.corrGetTexts()$references$Sellke_etal_2001) - } - if(options$effectSize){ mainTable$addColumnInfo(name = paste0(test, "_effect.size"), title = gettext("Effect size (Fisher's z)"), type = "number", overtitle = overtitle) mainTable$addColumnInfo(name = paste0(test, "_se.effect.size"), title = gettext("SE Effect size"), type = "number", overtitle = overtitle) } } } + + if(options$covariance){ + mainTable$addColumnInfo(name = "covariance", title = gettext("Covariance"), type = "number") + } } .corrTitlerer <- function(test, nTests){ @@ -252,13 +256,19 @@ CorrelationInternal <- function(jaspResults, dataset, options){ overtitle <- paste(vi, variables[vi], sep = ". ") if(options$sampleSize) { - mainTable$addColumnInfo(name = paste(.v(variables[vi]), "sample.size", sep = "_"), title = "n", + mainTable$addColumnInfo(name = paste(variables[vi], "sample.size", sep = "_"), title = "n", type = "integer", overtitle = overtitle) } for(ti in seq_along(tests)){ .corrInitCorrelationTableRowAsColumn(mainTable, options, variables[vi], testsTitles[ti], tests[ti], overtitle) } + + if(options$covariance) { + mainTable$addColumnInfo(name = paste(variables[vi], "covariance", sep = "_"), gettextf("Covariance"), + type = "number", overtitle = overtitle) + } + mainTable$setRowName(vi, .v(variables[vi])) } } @@ -290,11 +300,11 @@ CorrelationInternal <- function(jaspResults, dataset, options){ title = gettextf("Upper %s%% CI", 100*options$ciLevel), type = "number", overtitle = overtitle) } - + if(options$effectSize){ mainTable$addColumnInfo(name = sprintf(name, "effect.size"), title = gettextf("Effect size (Fisher's z)"), type = "number", overtitle = overtitle) mainTable$addColumnInfo(name = sprintf(name, "se.effect.size"), title = gettext("SE Effect size"), type = "number", overtitle = overtitle) - } + } } ### Compute results ---- @@ -343,7 +353,9 @@ CorrelationInternal <- function(jaspResults, dataset, options){ currentResults <- list() testErrors <- list() currentResults[['sample.size']] <- nrow(data) - + if (isFALSE(errors)) + currentResults[['covariance']] <- cov(x = data[,1], y = data[,2]) + # even if we do not want the specific tests results # we still want the output as NaN - to fill the jaspTables correctly # so we still loop over all tests - .corr.test() returns empty lists if isFALSE(compute) @@ -466,7 +478,6 @@ CorrelationInternal <- function(jaspResults, dataset, options){ statsNames <- c("estimate", "p.value", "lower.ci", "upper.ci", "vsmpr", "effect.size", "se.effect.size") alternative <- match.arg(alternative) - if(isFALSE(compute)){ result <- rep(NaN, length(statsNames)) names(result) <- statsNames @@ -513,9 +524,6 @@ CorrelationInternal <- function(jaspResults, dataset, options){ result$se.effect.size <- sqrt((2 / (n * (n - 1))) * (1 - (4 * (s1^2 / pi^2)) + (2 * (n - 2) * ((1/9) - (4 * (s2^2 / pi^2)))))) } - - - result <- unlist(result[stats], use.names = FALSE) names(result) <- statsNames } @@ -705,6 +713,7 @@ CorrelationInternal <- function(jaspResults, dataset, options){ testErrors <- lapply(corrResults, function(x) x[['testErrors']]) mainTable[['sample.size']] <- sapply(results, function(x) x[['sample.size']]) + mainTable[['covariance']] <- sapply(results, function(x) x[['covariance']]) # would be nice to be able to fill table cell-wise, i.e., mainTable[[row, col]] <- value colNames <- character() # this is for error footnotes diff --git a/R/generalizedlinearmodel.R b/R/generalizedlinearmodel.R index 161d35b3..725225fc 100644 --- a/R/generalizedlinearmodel.R +++ b/R/generalizedlinearmodel.R @@ -470,7 +470,14 @@ GeneralizedLinearModelInternal <- function(jaspResults, dataset = NULL, options, .glmOutlierTable(jaspResults, dataset, options, ready, position = 8, residType = "standardized deviance") .glmOutlierTable(jaspResults, dataset, options, ready, position = 8, residType = "studentized deviance") - .glmInfluenceTable(jaspResults, dataset, options, ready, position = 9) + .glmInfluenceTable(jaspResults[["diagnosticsContainer"]], + jaspResults[["glmModels"]][["object"]][["fullModel"]], + dataset, options, ready, position = 9) + + .regressionExportResiduals(jaspResults, + jaspResults[["glmModels"]][["object"]][["fullModel"]], + dataset, options, ready) + .glmMulticolliTable(jaspResults, dataset, options, ready, position = 10) return() @@ -639,53 +646,6 @@ GeneralizedLinearModelInternal <- function(jaspResults, dataset = NULL, options, } - -# Plots: Residuals Q-Q -.glmPlotResQQ <- function(jaspResults, dataset, options, ready, position) { - - plotNames <- c("devianceResidualQqPlot", "pearsonResidualQqPlot", "quantileResidualQqPlot") - if (!ready || !any(unlist(options[plotNames]))) - return() - - residNames <- c("deviance", "Pearson", "quantile") - - glmPlotResQQContainer <- createJaspContainer(gettext("Normal Q-Q Plots: Standardized Residuals")) - glmPlotResQQContainer$dependOn(optionsFromObject = jaspResults[["modelSummary"]], - options = c(plotNames, "seed", "setSeed")) - glmPlotResQQContainer$position <- position - jaspResults[["diagnosticsContainer"]][["glmPlotResQQ"]] <- glmPlotResQQContainer - - - if (!is.null(jaspResults[["glmModels"]])) { - glmFullModel <- jaspResults[["glmModels"]][["object"]][["fullModel"]] - for (i in 1:length(plotNames)) { - if (options[[plotNames[[i]]]]) { - .glmCreatePlotPlaceholder(glmPlotResQQContainer, - index = plotNames[[i]], - title = gettextf("Normal Q-Q plot: Standardized %1s residuals", residNames[[i]])) - - .glmInsertPlot(glmPlotResQQContainer[[plotNames[[i]]]], - .glmFillPlotResQQ, - residType = residNames[[i]], - model = glmFullModel, - family = options[["family"]]) - } - } - } - return() -} - -.glmFillPlotResQQ <- function(residType, model, family) { - - # compute residuals - stdResid <- .glmStdResidCompute(model = model, residType = residType, options = options) - - thePlot <- jaspGraphs::plotQQnorm(stdResid, ablineColor = "blue") - - return(thePlot) -} - - # Plots: Partial residuals .glmPlotResPartial <- function(jaspResults, dataset, options, ready, position) { if (!ready) @@ -894,112 +854,6 @@ GeneralizedLinearModelInternal <- function(jaspResults, dataset = NULL, options, } } - -# Table: Influential cases -.glmInfluenceTable <- function(jaspResults, dataset, options, ready, position) { - - tableOptionsOn <- c(options[["dfbetas"]], - options[["dffits"]], - options[["covarianceRatio"]], - options[["cooksDistance"]], - options[["leverage"]]) - - - if (!ready | !any(tableOptionsOn)) - return() - - - tableOptions <- c("dfbetas", "dffits", "covarianceRatio", "cooksDistance", "leverage") - tableOptionsClicked <- tableOptions[tableOptionsOn] - - if (is.null(jaspResults[["diagnosticsContainer"]][["influenceTable"]])) { - influenceTable <- createJaspTable(gettext("Table: Influential Cases")) - influenceTable$dependOn(optionsFromObject = jaspResults[["modelSummary"]], - options = tableOptions) - influenceTable$position <- position - influenceTable$showSpecifiedColumnsOnly <- TRUE - jaspResults[["diagnosticsContainer"]][["influenceTable"]] <- influenceTable - } - - tableOptionToColName <- function(x) { - switch(x, - "dfbetas" = "DFBETAS", - "dffits" = "DFFITS", - "covarianceRatio" = "Covariance Ratio", - "cooksDistance" = "Cook's Distance", - "leverage" = "Leverage") - } - - if (is.null(jaspResults[["glmModels"]])) { - for (option in tableOptionsClicked) { - colTitle <- tableOptionToColName(option) - jaspResults[["influenceTable"]]$addColumnInfo(name = option, title = gettext(colTitle), type = "number") - } - } else { - glmFullModel <- jaspResults[["glmModels"]][["object"]][["fullModel"]] - colNameList <- c() - jaspResults[["diagnosticsContainer"]][["influenceTable"]]$addColumnInfo(name = "caseN", title = "Case Number", type = "integer") - for (option in tableOptionsClicked) { - if (option == "dfbetas") { - predictors <- names(glmFullModel$coefficients) - for (predictor in predictors) { - dfbetasName <- gettextf("DFBETAS_%1s", predictor) - colNameList <- c(colNameList, dfbetasName) - if (predictor == "(Intercept)") - dfbetasTitle <- gettext("DFBETAS:Intercept") - else - dfbetasTitle <- gettextf("DFBETAS:%1s", gsub(":", "*", predictor)) - jaspResults[["diagnosticsContainer"]][["influenceTable"]]$addColumnInfo(name = dfbetasName, title = dfbetasTitle, type = "number") - } - } else { - colNameList <- c(colNameList, option) - colTitle <- tableOptionToColName(option) - jaspResults[["diagnosticsContainer"]][["influenceTable"]]$addColumnInfo(name = option, title = gettext(colTitle), type = "number") - } - } - .glmInfluenceTableFill(jaspResults, dataset, options, ready, model = glmFullModel, influenceMeasures = tableOptionsClicked, colNames = colNameList) - } -} - -.glmInfluenceTableFill <- function(jaspResults, dataset, options, ready, model, influenceMeasures, colNames) { - influenceRes <- influence.measures(model) - nDFBETAS <- length(names(model$coefficients)) - - optionToColInd <- function(x, nDFBETAS) { - switch(x, - "dfbetas" = 1:nDFBETAS, - "dffits" = (nDFBETAS+1), - "covarianceRatio" = (nDFBETAS+2), - "cooksDistance" = (nDFBETAS+3), - "leverage" = (nDFBETAS+4))} - - colInd <- c() - for (measure in influenceMeasures) { - colInd <- c(colInd, optionToColInd(measure, nDFBETAS)) - } - - influenceResData <- as.data.frame(influenceRes[["infmat"]][, colInd]) - names(influenceResData) <- colNames - caseN <- seq.int(nrow(influenceResData)) - influenceResData <- cbind(caseN, influenceResData) - - influenceResSig <- influenceRes[["is.inf"]][, colInd] - - if (length(colInd) > 1) { - influenceResDataFinal <- influenceResData[rowSums(influenceResSig) > 0, , drop = FALSE] - } else { - influenceResDataFinal <- influenceResData[influenceResSig > 0, , drop = FALSE] - } - - nRowInfluential <- nrow(influenceResDataFinal) - - if (nRowInfluential == 0) - jaspResults[["diagnosticsContainer"]][["influenceTable"]]$addFootnote(gettext("No influential cases found.")) - else { - jaspResults[["diagnosticsContainer"]][["influenceTable"]]$setData(influenceResDataFinal) - } -} - # Table: Multicollinearity .glmMulticolliTable <- function(jaspResults, dataset, options, ready, position) { diff --git a/R/glmCommonFunctions.R b/R/glmCommonFunctions.R deleted file mode 100644 index 6525f28b..00000000 --- a/R/glmCommonFunctions.R +++ /dev/null @@ -1,286 +0,0 @@ -# as explained in ?is.integer -.is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol - -# create GLM formula -.createGLMFormula <- function(options, nullModel = FALSE) { - - leftTerm <- options$dependent - modelTerms <- options$modelTerms - includeIntercept <- options$interceptTerm - offsetTerm <- options$offset - - if (includeIntercept) - rightTerms <- "1" - else - rightTerms <- "-1" - - if (offsetTerm != "") - rightTerms <- c(rightTerms, paste("offset(", offsetTerm, ")", sep = "")) - - if (length(modelTerms) == 0) { - f <- formula(paste(leftTerm, "~", rightTerms)) - } else { - - if (nullModel) { - for (i in seq_along(modelTerms)) { - nuisance <- modelTerms[[i]][["isNuisance"]] - if (!is.null(nuisance) && nuisance) { - term <- modelTerms[[i]][["components"]] - if (length(term) == 1) - rightTerms <- c(rightTerms, term) - else - rightTerms <- c(rightTerms, paste(term, collapse = ":")) - } - } - - } else { - for (i in seq_along(modelTerms)) { - term <- modelTerms[[i]][["components"]] - if (length(term) == 1) - rightTerms <- c(rightTerms, term) - else - rightTerms <- c(rightTerms, paste(term, collapse = ":")) - } - } - f <- formula(paste(leftTerm, "~", paste(rightTerms, collapse = "+"))) - } - return(f) -} - - -.glmComputeModel <- function(jaspResults, dataset, options) { - if (!is.null(jaspResults[["glmModels"]])) - return(jaspResults[["glmModels"]]$object) - - # make formulas for the full model and the null model - ff <- .createGLMFormula(options, nullModel = FALSE) - environment(ff) <- environment() - nf <- .createGLMFormula(options, nullModel = TRUE) - environment(nf) <- environment() - - weights <- dataset[[options[["weights"]]]] - - if (options$family != "other") { - # specify family and link - if (options$family == "bernoulli") { - family <- "binomial" - } - else if (options$family == "gamma") { - family <- "Gamma" - } - else if (options$family == "inverseGaussian") { - family <- "inverse.gaussian" - } - else { - family <- options$family - } - - familyLink <- eval(call(family, link = options$link)) - # compute full and null models - fullModel <- try(stats::glm(ff, - family = familyLink, - data = dataset, - weights = weights)) - nullModel <- try(stats::glm(nf, - family = familyLink, - data = dataset, - weights = weights)) - } else { - - if (options$otherGlmModel == "multinomialLogistic") { - fullModel <- try(VGAM::vglm(ff, - family = VGAM::multinomial(), - data = dataset, - weights = weights)) - nullModel <- try(VGAM::vglm(nf, - family = VGAM::multinomial(), - data = dataset, - weights = weights)) - } - - if (options$otherGlmModel == "ordinalLogistic") { - fullModel <- try(VGAM::vglm(ff, - family = VGAM::cumulative(link = "logitlink", parallel = TRUE), - data = dataset, - weights = weights)) - nullModel <- try(VGAM::vglm(nf, - family = VGAM::cumulative(link = "logitlink", parallel = TRUE), - data = dataset, - weights = weights)) - } - - if (options$otherGlmModel == "firthLogistic") { - fullModel <- try(logistf::logistf(ff, - data = dataset, - pl = TRUE, - firth = TRUE, - weights = weights)) - nullModel <- try(logistf::logistf(nf, - data = dataset, - pl = TRUE, - firth = TRUE, - weights = weights)) - } - } - - if (jaspBase::isTryError(fullModel)) { - msg <- .glmFitErrorMessageHelper(fullModel) - msg <- gettextf("The full model could not be fitted to the data, with the following error message: '%s'", msg) - jaspBase::.quitAnalysis(msg) - } - - if (jaspBase::isTryError(nullModel)) { - msg <- .glmFitErrorMessageHelper(nullModel) - msg <- gettextf("The null model could not be fitted to the data, with the following error message: '%s'", msg) - jaspBase::.quitAnalysis(msg) - } - - glmModels <- list("nullModel" = nullModel, - "fullModel" = fullModel) - # combine both models - jaspResults[["glmModels"]] <- createJaspState() - jaspResults[["glmModels"]]$dependOn(optionsFromObject = jaspResults[["modelSummary"]]) - jaspResults[["glmModels"]]$object <- glmModels - - return(glmModels) -} - -.glmFitErrorMessageHelper <- function(tryError) { - msg <- jaspBase::.extractErrorMessage(tryError) - if(msg == "cannot find valid starting values: please specify some") - msg <- gettext("Could not find valid starting values. Check feasibility of the model to fit the data.") - - return(msg) -} - - -.hasNuisance <- function(options) { - return(any(sapply(options$modelTerms, function(x) x[["isNuisance"]]))) -} - -.glmCreatePlotPlaceholder <- function(container, index, title) { - jaspPlot <- createJaspPlot(title = title) - jaspPlot$status <- "running" - container[[index]] <- jaspPlot - return() -} - - -.glmStdResidCompute <- function(model, residType, options) { - if (residType == "deviance") { - stdResid <- rstandard(model) - } - else if (residType == "Pearson") { - phiEst <- summary(model)[["dispersion"]] - resid <- resid(model, type="pearson") - stdResid <- resid / sqrt(phiEst * (1 - hatvalues(model))) - } - else if (residType == "quantile") { - jaspBase::.setSeedJASP(options) - resid <- statmod::qresid(model) - stdResid <- resid / sqrt(1 - hatvalues(model)) - } else { - stdResid <- NULL - } - - return(stdResid) -} - -.glmInsertPlot <- function(jaspPlot, func, ...) { - p <- try(func(...)) - - if (inherits(p, "try-error")) { - errorMessage <- .extractErrorMessage(p) - jaspPlot$setError(gettextf("Plotting is not possible: %s", errorMessage)) - } else { - jaspPlot$plotObject <- p - jaspPlot$status <- "complete" - } -} - -.constInfoTransform <- function(family, x) { - switch(family, - "bernoulli" = 1/(sin(sqrt(x))), - "binomial" = 1/(sin(sqrt(x))), - "poisson" = sqrt(x), - "Gamma" = log(x), - "inverse.gaussian" = 1/sqrt(x), - "gaussian" = x) -} - -.constInfoTransName <- function(family) { - switch(family, - "bernoulli" = expression(sin^-1 * sqrt("Fitted values")), - "binomial" = expression(sin^-1 * sqrt("Fitted values")), - "poisson" = expression(sqrt("Fitted values")), - "Gamma" = expression(log("Fitted values")), - "inverse.gaussian" = expression(1/sqrt("Fitted Values")), - "gaussian" = "Fitted values") -} - -.capitalize <- function(x) { - x_new <- paste(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)), sep="") - return(x_new) -} - -# function for multicollineary statistics, taken from the source code of car::vif -# car version: 3.0-13 see https://rdrr.io/cran/car/src/R/vif.R -# the reason is to remove the dependency on the car package. When importing the entire car package, -# JASP crashes when loading the regression module (at least on Windows 10). - -.vif.default <- function(mod, ...) { - # modified to fix bug: https://github.com/jasp-stats/jasp-test-release/issues/2487 - if (length(coef(mod)) == 0) { - stop(gettext("There are no predictors to test for multicollinearity")) - } - # end modification - if (any(is.na(coef(mod)))) - stop ("there are aliased coefficients in the model") - v <- vcov(mod) - assign <- attr(model.matrix(mod), "assign") - if (names(coefficients(mod)[1]) == "(Intercept)") { - v <- v[-1, -1] - assign <- assign[-1] - } - else warning("No intercept: vifs may not be sensible.") - 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 <- matrix(0, n.terms, 3) - rownames(result) <- terms - colnames(result) <- c("GVIF", "Df", "GVIF^(1/(2*Df))") - for (term in 1:n.terms) { - subs <- which(assign == term) - result[term, 1] <- det(as.matrix(R[subs, subs])) * - det(as.matrix(R[-subs, -subs])) / detR - result[term, 2] <- length(subs) - } - if (all(result[, 2] == 1)) result <- result[, 1] - else result[, 3] <- result[, 1]^(1/(2 * result[, 2])) - result -} - -# message for estimated marginal means table -.emmMessageTestNull <- function(value) gettextf("P-values correspond to test of null hypothesis against %s.", value) -.emmMessageAveragedOver <- function(terms) gettextf("Results are averaged over the levels of: %s.",paste(terms, collapse = ", ")) -.messagePvalAdjustment <- function(adjustment) { - if (adjustment == "none") { - return(gettext("P-values are not adjusted.")) - } - adjustment <- switch(adjustment, - "holm" = gettext("Holm"), - "hommel" = gettext("Homel"), - "hochberg" = gettext("Hochberg"), - "mvt" = gettext("Multivariate-t"), - "tukey" = gettext("Tukey"), - "BH" = gettext("Benjamini-Hochberg"), - "BY" = gettext("Benjamini-Yekutieli"), - "scheffe" = gettext("Scheffé"), - "sidak" = gettext("Sidak"), - "dunnettx" = gettext("Dunnett"), - "bonferroni" = gettext("Bonferroni") - ) - return(gettextf("P-values are adjusted using %s adjustment.", adjustment)) -} diff --git a/R/regressionlinear.R b/R/regressionlinear.R index f33ac7f8..f7e8ae41 100755 --- a/R/regressionlinear.R +++ b/R/regressionlinear.R @@ -51,9 +51,11 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { # these output elements show statistics of the "final model" (lm fit with all predictors in enter method and last lm fit in stepping methods) finalModel <- model[[length(model)]] - if (options$residualCasewiseDiagnostic && is.null(modelContainer[["casewiseTable"]])) - .linregCreateCasewiseDiagnosticsTable(modelContainer, finalModel, options, position = 9) - + if (options$residualCasewiseDiagnostic && is.null(modelContainer[["influenceTable"]])) + .glmInfluenceTable(modelContainer, finalModel$fit, dataset, options, ready = ready, position = 9, linRegAnalysis = TRUE) + .regressionExportResiduals(modelContainer, finalModel$fit, dataset, options, ready = ready) + + if (options$residualStatistic && is.null(modelContainer[["residualsTable"]])) .linregCreateResidualsTable(modelContainer, finalModel, options, position = 10) @@ -108,12 +110,12 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { .linregCheckIfFactorWithMoreLevels <- function(var) { # Custom function to check if a variable is a factor with more than 2 levels - is.factor(var) && nlevels(var) > 2 + return(is.factor(var) && nlevels(var) > 2) } .linregCheckIfInteractionWithFactors <- function(modelTerm, factorVariables) { # Custom function to check if interaction contains more than 1 factor - sum(modelTerm[["components"]] %in% factorVariables) > 1 + return(sum(modelTerm[["components"]] %in% factorVariables) > 1) } .linregCheckErrors <- function(dataset, options) { @@ -147,16 +149,13 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { } defaultTarget <- c(options$dependent, unlist(options$covariates)) - .hasErrors(dataset, type = c("infinity", "variance", "observations", "modelInteractions", "varCovData"), + .hasErrors(dataset, type = c("infinity", "variance", "observations", "varCovData"), custom = stepwiseProcedureChecks, custom.target = defaultTarget, observations.amount = "< 2", observations.target = defaultTarget, - modelInteractions.modelTerms = options$modelTerms, - modelInteractions.target = defaultTarget, - varCovData.target = unlist(options$covariates), varCovData.corFun = stats::cov, @@ -191,7 +190,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { else summaryTable <- createJaspTable(gettextf("Model Summary - %s", options[['dependent']])) - summaryTable$dependOn(c("residualDurbinWatson", "rSquaredChange")) + summaryTable$dependOn(c("residualDurbinWatson", "rSquaredChange", "fChange", "modelAICBIC")) summaryTable$position <- position summaryTable$showSpecifiedColumnsOnly <- TRUE @@ -201,9 +200,16 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { summaryTable$addColumnInfo(name = "adjR2", title = gettextf("Adjusted R%s", "\u00B2"), type = "number", format = "dp:3") summaryTable$addColumnInfo(name = "RMSE", title = gettext("RMSE"), type = "number") - if (options$rSquaredChange) { - summaryTable$addColumnInfo(name = "R2c", title = gettextf("R%s Change", "\u00B2"), type = "number", format = "dp:3") - summaryTable$addColumnInfo(name = "Fc", title = gettext("F Change"), type = "number") + if (options$modelAICBIC) { + summaryTable$addColumnInfo(name = "AIC", title = gettext("AIC"), type = "number", format = "dp:3") + summaryTable$addColumnInfo(name = "BIC", title = gettext("BIC"), type = "number", format = "dp:3") + } + + if (options$rSquaredChange || options$fChange) { + if (options$rSquaredChange) + summaryTable$addColumnInfo(name = "R2c", title = gettextf("R%s Change", "\u00B2"), type = "number", format = "dp:3") + if (options$fChange) + summaryTable$addColumnInfo(name = "Fc", title = gettext("F Change"), type = "number") summaryTable$addColumnInfo(name = "df1", title = gettext("df1"), type = "integer") summaryTable$addColumnInfo(name = "df2", title = gettext("df2"), type = "integer") summaryTable$addColumnInfo(name = "p", title = gettext("p"), type = "pvalue") @@ -219,7 +225,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { summaryTable$addFootnote(message = gettext("p-value for Durbin-Watson test is unavailable for weighted regression.")) } - .linregAddPredictorsInNullFootnote(summaryTable, options$modelTerms) + .linregAddPredictorsInNullFootnote(summaryTable, options$modelTerms[[1]][["components"]]) if (!is.null(model)) { if (length(model) == 1 && length(model[[1]]$predictors) == 0 && !options$interceptTerm) @@ -242,6 +248,8 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { R = as.numeric(sqrt(lmSummary$r.squared)), R2 = as.numeric(lmSummary$r.squared), adjR2 = as.numeric(lmSummary$adj.r.squared), + AIC = as.numeric(AIC(model[[i]][["fit"]])), + BIC = as.numeric(BIC(model[[i]][["fit"]])), RMSE = as.numeric(lmSummary$sigma), R2c = rSquareChange$R2c, Fc = rSquareChange$Fc, @@ -269,7 +277,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { anovaTable$addColumnInfo(name = "F", title = gettext("F"), type = "number") anovaTable$addColumnInfo(name = "p", title = gettext("p"), type = "pvalue") - .linregAddPredictorsInNullFootnote(anovaTable, options$modelTerms) + .linregAddPredictorsInNullFootnote(anovaTable, options$modelTerms[[1]][["components"]]) .linregAddVovkSellke(anovaTable, options$vovkSellke) if (!is.null(model)) { @@ -299,7 +307,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { coeffTable <- createJaspTable(gettext("Coefficients")) coeffTable$dependOn(c("coefficientEstimate", "coefficientCi", "coefficientCiLevel", - "collinearityDiagnostic", "vovkSellke")) + "collinearityStatistic", "vovkSellke")) coeffTable$position <- position coeffTable$showSpecifiedColumnsOnly <- TRUE @@ -319,7 +327,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { coeffTable$addColumnInfo(name = "upper", title = gettext("Upper"), type = "number", overtitle = overtitle) } - if (options$collinearityDiagnostic) { + if (options$collinearityStatistic) { overtitle <- gettext("Collinearity Statistics") coeffTable$addColumnInfo(name = "tolerance", title = gettext("Tolerance"), type = "number", format = "dp:3", overtitle = overtitle) coeffTable$addColumnInfo(name = "VIF", title = gettext("VIF"), type = "number", overtitle = overtitle) @@ -579,37 +587,6 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { } } -.linregCreateCasewiseDiagnosticsTable <- function(modelContainer, finalModel, options, position) { - caseDiagTable <- createJaspTable(gettext("Casewise Diagnostics")) - caseDiagTable$dependOn(c("residualCasewiseDiagnostic", "residualCasewiseDiagnosticType", - "residualCasewiseDiagnosticZThreshold", "residualCasewiseDiagnosticCooksDistanceThreshold")) - caseDiagTable$position <- position - - caseDiagTable$addColumnInfo(name = "caseNumber", title = gettext("Case Number"), type = "integer") - caseDiagTable$addColumnInfo(name = "stdResidual", title = gettext("Std. Residual"), type = "number", format = "dp:3") - caseDiagTable$addColumnInfo(name = "dependent", title = options$dependent, type = "number") - caseDiagTable$addColumnInfo(name = "predicted", title = gettext("Predicted Value"), type = "number") - caseDiagTable$addColumnInfo(name = "residual", title = gettext("Residual"), type = "number") - caseDiagTable$addColumnInfo(name = "cooksD", title = gettext("Cook's Distance"), type = "number", format = "dp:3") - - if (!is.null(finalModel)) { - caseDiagData <- .linregGetCasewiseDiagnostics(finalModel$fit, options) - caseDiagTable$setData(caseDiagData) - - if (length(caseDiagData) == 0) { - message <- switch( - options[["residualCasewiseDiagnosticType"]], - cooksDistance = gettextf("No cases where |Cook's distance| > %s", options[["residualCasewiseDiagnosticCooksDistanceThreshold"]]), - outliersOutside = gettextf("No cases where |Standard residual| > %s", options[["residualCasewiseDiagnosticZThreshold"]]), - gettextf("No cases to show") - ) - caseDiagTable$addFootnote(message = message) - } - } - - modelContainer[["casewiseTable"]] <- caseDiagTable -} - .linregCreateResidualsTable <- function(modelContainer, finalModel, options, position) { residualsTable <- createJaspTable(gettext("Residuals Statistics")) residualsTable$dependOn("residualStatistic") @@ -724,12 +701,14 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { if (!is.null(finalModel) && !is.null(finalModel$fit)) { fit <- finalModel$fit - .linregInsertPlot(residQQPlot, .linregPlotQQresiduals, res = residuals(fit) / sd(residuals(fit))) + + .linregInsertPlot(residQQPlot, .glmFillPlotResQQ, model = fit, + residType = "deviance", options = options) } } .linregCreatePartialPlots <- function(modelContainer, dataset, options, position) { - predictors <- .linregGetPredictors(options$modelTerms, encoded = TRUE) + predictors <- .linregGetPredictors(options$modelTerms) title <- ngettext(length(predictors), "Partial Regression Plot", "Partial Regression Plots") @@ -739,7 +718,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { partialPlotContainer$position <- position modelContainer[["partialPlotContainer"]] <- partialPlotContainer - predictors <- .linregGetPredictors(options$modelTerms, encoded = TRUE) + predictors <- .linregGetPredictors(options$modelTerms[[length(options$modelTerms)]][["components"]]) if (any(.linregIsInteraction(predictors))) { .linregCreatePlotPlaceholder(partialPlotContainer, index = "placeholder", title = "") partialPlotContainer$setError(gettext("Partial plots are not supported for models containing interaction terms")) @@ -811,24 +790,26 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { model <- .linregCalcDurBinWatsonTestResults(modelContainer, model, options) return(model) } - + nModels <- length(options$modelTerms) dependent <- options$dependent - predictorsInNull <- .linregGetPredictors(options$modelTerms, modelType = "null") - predictors <- .linregGetPredictors(options$modelTerms, modelType = "alternative") # these include the null terms + predictorsInNull <- .linregGetPredictors(options$modelTerms[[1]][["components"]], modelType = "null") + predictorsInFull <- .linregGetPredictors(options$modelTerms[[nModels]][["components"]], modelType = "alternative") # these include the null terms + + if (options$weights != "") weights <- dataset[[options$weights]] else weights <- rep(1, length(dataset[[dependent]])) - - if (options$method %in% c("backward", "forward", "stepwise") && length(predictors) > 0) - model <- .linregGetModelSteppingMethod(dependent, predictors, predictorsInNull, dataset, options, weights) + + if (options$method %in% c("backward", "forward", "stepwise") && length(predictorsInFull) > 0) + model <- .linregGetModelSteppingMethod(dependent, predictorsInFull, predictorsInNull, dataset, options, weights) else - model <- .linregGetModelEnterMethod(dependent, predictors, predictorsInNull, dataset, options, weights) + model <- .linregGetModelEnterMethod(dependent, modelTerms = options[["modelTerms"]], dataset, options, weights) for (i in seq_along(model)) { singleModel <- model[[i]] - model[[i]][["title"]] <- .linregGetModelTitle(singleModel$predictors, predictorsInNull, options$method, i) + model[[i]][["title"]] <- gettextf("M%s", intToUtf8(0x2080 + i - 1, multiple = FALSE)) # singleModel[["title"]] model[[i]][["summary"]] <- .linregGetSummary(singleModel$fit) model[[i]][["rSquareChange"]] <- .linregGetrSquaredChange(singleModel$fit, i, model[1:i], options) } @@ -868,27 +849,15 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { return(model) } -.linregGetModelEnterMethod <- function(dependent, predictors, predictorsInNull, dataset, options, weights) { +.linregGetModelEnterMethod <- function(dependent, modelTerms, dataset, options, weights) { model <- list() - formulaNull <- NULL - if (length(predictorsInNull) > 0) - formulaNull <- .linregGetFormula(dependent, predictorsInNull, options$interceptTerm) - else if (options$interceptTerm) - formulaNull <- .linregGetFormula(dependent, NULL, TRUE) - - formula <- NULL - if (length(predictors) > 0 && !all(predictors %in% predictorsInNull)) - formula <- .linregGetFormula(dependent, predictors, options$interceptTerm) - - if (!is.null(formulaNull)) { - fitNull <- stats::lm(formulaNull, data = dataset, weights = weights, x = TRUE) - model[[1]] <- list(fit = fitNull, predictors = predictorsInNull, title = gettextf("H%s", "\u2080")) - } - - if (!is.null(formula)) { + for (thisModel in modelTerms) { + thisModelTerms <- .linregGetPredictors(thisModel[["components"]]) + formula <- .linregGetFormula(dependent, thisModelTerms, options$interceptTerm) fit <- stats::lm(formula, data = dataset, weights = weights, x = TRUE) - model[[length(model) + 1]] <- list(fit = fit, predictors = predictors, title = gettextf("H%s", "\u2081")) + model[[length(model) + 1]] <- list(fit = fit, predictors = thisModelTerms, title = gettext(thisModel[["title"]])) + } return(model) @@ -1504,6 +1473,8 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { residualsAll <- residuals(fit) stdPredictedValuesAll <- (predictedValuesAll - mean(predictedValuesAll)) / sd(predictedValuesAll) stdResidualsAll <- rstandard(fit) + # stdResidualsAll <- statmod::qresid(fit) + # stdResidualsAll <- rstudent(fit cooksDAll <- cooks.distance(fit) if (options$residualCasewiseDiagnosticType == "cooksDistance") @@ -1755,65 +1726,18 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { return(p) } -.linregPlotQQresiduals <- function(res = NULL, xlab = gettext("Theoretical Quantiles"), ylab= gettext("Standardized Residuals"), cexPoints= 1.3, cexXAxis= 1.3, cexYAxis= 1.3, lwd= 2, lwdAxis=1.2) { - d <- data.frame(qqnorm(res, plot.it = FALSE)) - d <- na.omit(d) - xVar <- d$x - yVar <- d$y - - xlow <- min(pretty(xVar)) - xhigh <- max(pretty(xVar)) - xticks <- pretty(c(xlow, xhigh)) - - ylow <- min(pretty(yVar)) - yhigh <- max(pretty(yVar)) - yticks <- pretty(c(ylow, yhigh)) - - yLabs <- vector("character", length(yticks)) - - for (i in seq_along(yticks)) { - if (yticks[i] < 10^6) - yLabs[i] <- format(yticks[i], digits= 3, scientific = FALSE) - else - yLabs[i] <- format(yticks[i], digits= 3, scientific = TRUE) - } - - p <- ggplot2::ggplot() + - ggplot2::scale_x_continuous(name = gettext("Theoretical Quantiles"), breaks = xticks) + - ggplot2::scale_y_continuous(name = gettext("Standardized Residuals"), breaks = xticks) +# TODO: why are xticks used here even though yticks exists? - ggplot2::geom_line( - data = data.frame(x = c(min(xticks), max(xticks)), y = c(min(xticks), max(xticks))), - mapping = ggplot2::aes(x = x, y = y), - col = "darkred", size = 1 - ) + - jaspGraphs::geom_point( - data = data.frame(x = xVar, y = yVar), - mapping = ggplot2::aes(x = x, y = y) - ) + - jaspGraphs::geom_rangeframe() + - jaspGraphs::themeJaspRaw(axis.title.cex = 1.2) - - return(p) -} - -.linregGetPredictors <- function(modelTerms, modelType = "alternative", encoded = TRUE) { +.linregGetPredictors <- function(modelTerms, modelType = "alternative") { if (!is.character(modelType) || !modelType %in% c("alternative", "null")) stop(gettext("Unknown value provided for modelType, possible values: `alternative`, `null`")) predictors <- NULL for (i in seq_along(modelTerms)) { - components <- unlist(modelTerms[[i]]$components) - if (encoded) - components <- .v(components) + components <- modelTerms[[i]] predictor <- paste0(components, collapse = ":") - if (modelType == "alternative") { - predictors <- c(predictors, predictor) - } else if (modelType == "null") { - isNuisance <- modelTerms[[i]]$isNuisance - if (isNuisance) - predictors <- c(predictors, predictor) - } + + predictors <- c(predictors, predictor) + } return(predictors) @@ -1833,7 +1757,7 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { if (is.null(predictors) && includeConstant == FALSE) stop(gettext("We need at least one predictor, or an intercept to make a formula")) - if (is.null(predictors)) + if (length(predictors) == 0) formula <- paste(dependent, "~", "1") else if (includeConstant) formula <- paste(dependent, "~", paste(predictors, collapse = "+")) @@ -1843,17 +1767,6 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { return(as.formula(formula, env = parent.frame(1))) } -.linregGetModelTitle <- function(predictors, predictorsInNull, method, index) { - modelTitle <- index - if (method == "enter") { - if (index == 1 && (length(predictors) == 0 || all(predictors %in% predictorsInNull))) - modelTitle <- gettextf("H%s", "\u2080") - else - modelTitle <- gettextf("H%s", "\u2081") - } - - return(modelTitle) -} .linregIsInteraction <- function(predictor) { grepl(":", predictor) @@ -1899,9 +1812,8 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { } .linregAddPredictorsInNullFootnote <- function(jaspTable, modelTerms) { - containsNuisance <- .reglinTermsContainNuisance(modelTerms) - if (containsNuisance) { - predictorsInNull <- .linregGetPredictors(modelTerms, modelType = "null", encoded = FALSE) + if (length(modelTerms) > 0) { + predictorsInNull <- .linregGetPredictors(modelTerms, modelType = "null") jaspTable$addFootnote(message = gettextf("Null model includes %s", paste(predictorsInNull, collapse = ", "), sep = "")) } } @@ -1926,10 +1838,11 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) { } .linregGetIndicesOfModelsWithPredictors <- function(model, options) { - predictorsInNull <- .linregGetPredictors(options$modelTerms, modelType = "null") + predictorsInNull <- model[[1]]$predictors indices <- seq_along(model) + if (options$method == "enter") { - if (length(model) >= 1 && options$interceptTerm && is.null(predictorsInNull)) + if (length(model) >= 1 && options$interceptTerm && length(predictorsInNull) == 0) indices <- indices[-1] } else { for (i in seq_along(model)) diff --git a/R/regressionlinearbayesian.R b/R/regressionlinearbayesian.R index f98c23cc..791aa27c 100644 --- a/R/regressionlinearbayesian.R +++ b/R/regressionlinearbayesian.R @@ -762,7 +762,7 @@ for sparse regression when there are more covariates than observations (Castillo p <- try({ x <- fitted(basregModel, estimator = "BMA") y <- basregModel$Y - x - jaspGraphs::plotQQnorm(y) + jaspGraphs::plotQQnorm(y, ablineColor = "darkred") }) if (isTryError(p)) { diff --git a/inst/help/Correlation.md b/inst/help/Correlation.md index 11f44767..c0e5d65d 100755 --- a/inst/help/Correlation.md +++ b/inst/help/Correlation.md @@ -43,6 +43,7 @@ The Correlation analysis allows estimation of the population correlation, as wel - Vovk-Selke maximum p-ratio: The bound 1/(-e p log(p)) is derived from the shape of the p-value distribution. Under the null hypothesis (H0) it is uniform (0,1), and under the alternative (H1) it is decreasing in p, e.g., a beta (α, 1) distribution, where 0 < α < 1. The Vovk-Sellke MPR is obtained by choosing the shape α of the distribution under H1 such that the obtained p-value is maximally diagnostic. The value is then the ratio of the densities at point p under H0 and H1. For example, if the two-sided p-value equals .05, the Vovk-Sellke MPR equals 2.46, indicating that this p-value is at most 2.46 times more likely to occur under H1 than under H0. - Effect size (Fisher's z): The Fisher transformed effect size with standard error. - Sample size: The number of complete observations for a given pair of variables. +- Covariance: The covariance between each pair of variables. #### Plots - Scatter plots: Display a scatter plots for each possible combination of the selected variables. In a matrix format, these are placed above the diagonal. diff --git a/inst/help/GeneralizedLinearModel.md b/inst/help/GeneralizedLinearModel.md index 41f108c2..b3b5235c 100644 --- a/inst/help/GeneralizedLinearModel.md +++ b/inst/help/GeneralizedLinearModel.md @@ -79,12 +79,18 @@ The following table summarized the available distributions (also called families - Top n standardized quantile residuals - Top n standardized deviance residuals - Top n studentized deviance residuals -- Show influential cases: A table showing influential observations (i.e. outliers with high leverage) according to the selected measure(s). The column "Case Number" refers to the row number of the observation in the data set. The cut-offs for determining whether a case is influential are listed as follows. - - DFBETAS: When the absolute value of DFBETAS is greater than 1. - - DFFITS: When the absolute value of DFFITS is greater than 3 * sqrt(k/(n-k)) where k refers to the number of parameters in the model and n refers to the sample size. - - Covariance ratio: When the covariance ratio is greater than 3 * k/(n-k). - - Cook's distance: When Cook's distance exceeds the 50th percentile of the F distribution with (k, n-k) degrees of freedom. - - Leverages: When the leverages are greater than 3 * k/n. +- Residuals: + - Casewise diagnostic: Casewise and summarized diagnostics for the residuals. + - Standard residual > 3: Outliers outside x standard deviations: Display diagnostics for cases where the absolute value of the standardized residual is larger than x; default is x=3. + - Cook's distance > 1 : Display diagnostics for cases where the value of Cook’s distance is larger than x; default is x = 1. + - All cases: Display diagnostics for all cases. + - Cases are marked as influential in the table, according to the following thresholds: + - DFBETAS: When the absolute value of DFBETAS is greater than 1. + - DFFITS: When the absolute value of DFFITS is greater than 3 * sqrt(k/(n-k)) where k refers to the number of parameters in the model and n refers to the sample size. + - Covariance ratio: When the covariance ratio is greater than 3 * k/(n-k). + - Cook's distance: When Cook's distance exceeds the 50th percentile of the F distribution with (k, n-k) degrees of freedom. + - Leverages: When the leverages are greater than 3 * k/n. + - Append residuals to data: Save the residuals of the most complex model as a new column in the data file. - Multicollinearity: A table showing multicollinearity diagnostics of the model. The choices of measures are as follows. - Tolerance - VIF: Variance Inflation Factor. diff --git a/inst/help/RegressionLinear.md b/inst/help/RegressionLinear.md index 1e68b090..442b7fa2 100755 --- a/inst/help/RegressionLinear.md +++ b/inst/help/RegressionLinear.md @@ -39,22 +39,38 @@ Linear regression allows the user to model a linear relationship between one or - Estimates: Unstandardized and standardized coefficient estimates, standard errors, t-values, and their corresponding p-values. - From `...` bootstraps: By selecting this option, bootstrapped estimation is applied. By default, the number of replications is set to 1000. This can be changed into the desired number. - Confidence Intervals: By selecting this option, confidence intervals for the estimated mean difference will be included. By default the confidence level is set to 95%. This can be changed into the desired percentage. - - Covariance matrix: Display the covariance matrix of the predictor variables, per model. + - Tolerance and VIF: Display Tolerance and Variance Inflation Factor for each predictor in the model to assess multicollinearity. - Vovk-Sellke Maximum *p*-Ratio: The bound 1/(-e *p* log(*p*)) is derived from the shape of the *p*-value distribution. Under the null hypothesis (H0) it is uniform(0,1), and under the alternative (H1) it is decreasing in *p*, e.g., a beta(α, 1) distribution, where 0 < α < 1. The Vovk-Sellke MPR is obtained by choosing the shape α of the distribution under H1 such that the obtained *p*-value is *maximally diagnostic*. The value is then the ratio of the densities at point *p* under H0 and H1. For example, if the two-sided *p*-value equals .05, the Vovk-Sellke MPR equals 2.46, indicating that this *p*-value is at most 2.46 times more likely to occur under H1 than under H0. + +- Model Summary: + - R squared change: Change in R squared between the different steps in Backward, Forward, and Stepwise regression, with corresponding significance test (i.e., df1, df2, p-value). + - F change: Change in F between the different steps in Backward, Forward, and Stepwise regression, with corresponding significance test (i.e., df1, df2, p-value). + - AIC and BIC: Display Akaike Information Criterion and Bayesian Information Criterion. + - Durbin-Watson: Durbin-Watson statistic to test the autocorrelation of the residuals. + +- Display: - Model fit: Separate ANOVA table for each model (i.e., each step in Backward, Forward, and Stepwise regression). - - R squared change: Change in R squared between the different steps in Backward, Forward, and Stepwise regression, with corresponding significance test (i.e., F change value, df1, df2, p-value). - Descriptives: Samples size, sample mean, sample standard deviation, and standard error of the mean. - Part and partial correlations: Semipartial and partial correlations. + - Coefficient covariance matrix: Display the covariance matrix of the predictor variables, per model. - Collinearity diagnostics: Collinearity statistics, eigenvalues, condition indices, and variance proportions. - Residuals: - Statistics: Display descriptive statistics of the residuals and predicted values. - - Durbin-Watson: Durbin-Watson statistic to test the autocorrelation of the residuals. - Casewise diagnostic: Casewise and summarized diagnostics for the residuals. - Standard residual > 3: Outliers outside x standard deviations: Display diagnostics for cases where the absolute value of the standardized residual is larger than x; default is x=3. - Cook's distance > 1 : Display diagnostics for cases where the value of Cook’s distance is larger than x; default is x = 1. - All cases: Display diagnostics for all cases. + - Cases are marked as influential in the table, according to the following thresholds: + - DFBETAS: When the absolute value of DFBETAS is greater than 1. + - DFFITS: When the absolute value of DFFITS is greater than 3 * sqrt(k/(n-k)) where k refers to the number of parameters in the model and n refers to the sample size. + - Covariance ratio: When the covariance ratio is greater than 3 * k/(n-k). + - Cook's distance: When Cook's distance exceeds the 50th percentile of the F distribution with (k, n-k) degrees of freedom. + - Leverages: When the leverages are greater than 3 * k/n. + - Append residuals to data: Save the residuals of the most complex model as a new column in the data file. + + ### Method Specification diff --git a/inst/qml/Correlation.qml b/inst/qml/Correlation.qml index 3e449873..e07afdd7 100644 --- a/inst/qml/Correlation.qml +++ b/inst/qml/Correlation.qml @@ -68,6 +68,8 @@ Form CheckBox { name: "vovkSellke"; label: qsTr("Vovk-Sellke maximum p-ratio") } CheckBox { name: "effectSize"; label: qsTr("Effect size (Fisher's z)") } CheckBox { name: "sampleSize"; label: qsTr("Sample size") } + CheckBox { name: "covariance"; label: qsTr("Covariance") } + } diff --git a/inst/qml/GeneralizedLinearModel.qml b/inst/qml/GeneralizedLinearModel.qml index f21b2094..090d68ac 100644 --- a/inst/qml/GeneralizedLinearModel.qml +++ b/inst/qml/GeneralizedLinearModel.qml @@ -20,11 +20,15 @@ import QtQuick import QtQuick.Layouts import JASP import JASP.Controls -import "./common" as GLM +import "./common" as Common // All Analysis forms must be built with the From QML item Form { + id: form + property int analysis: Common.Type.Analysis.GLM + property int framework: Common.Type.Framework.Classical + Formula { lhs: "dependent" @@ -32,7 +36,7 @@ Form userMustSpecify: "covariates" } - GLM.GlmInputComponent { id: input} + Common.GlmInputComponent {} Section { @@ -94,7 +98,9 @@ Form title: qsTr("Diagnostics") enabled: input.otherFamilyNotSelected - GLM.GlmResidualAnalysisPlotsComponent {} + Common.OutlierComponent { id: outlierComponentt} + + Common.GlmResidualAnalysisPlotsComponent {} Group { @@ -102,7 +108,7 @@ Form CheckBox { name: "quantileResidualOutlierTable" - label: qsTr("Standardized quantile residuals: top") + label: qsTr("Standardized quantile residuals : top") childrenOnSameRow: true IntegerField { name: "quantileResidualOutlierTableTopN"; defaultValue: 3 } } @@ -122,15 +128,7 @@ Form } } - Group - { - title: qsTr("Show Influential Cases") - CheckBox { name: "dfbetas"; label: qsTr("DFBETAS") } - CheckBox { name: "dffits"; label: qsTr("DFFITS") } - CheckBox { name: "covarianceRatio"; label: qsTr("Covariance ratio") } - CheckBox { name: "cooksDistance"; label: qsTr("Cook's distance") } - CheckBox { name: "leverage"; label: qsTr("Leverages") } - } + Group { @@ -140,7 +138,7 @@ Form } } - GLM.EmmComponent { enabled: input.otherFamilyNotSelected } + Common.EmmComponent { enabled: input.otherFamilyNotSelected } Section { diff --git a/inst/qml/RegressionLinear.qml b/inst/qml/RegressionLinear.qml index 9f82fec6..d779f171 100644 --- a/inst/qml/RegressionLinear.qml +++ b/inst/qml/RegressionLinear.qml @@ -19,9 +19,14 @@ import QtQuick import QtQuick.Layouts import JASP import JASP.Controls +import "./common" as Common Form { + id: form + property int analysis: Common.Type.Analysis.LinearRegression + property int framework: Common.Type.Framework.Classical + Formula { lhs: "dependent" @@ -49,22 +54,31 @@ Form AssignedVariablesList { name: "weights"; title: qsTr("WLS Weights (optional)"); allowedColumns: ["scale"]; singleVariable: true } } - Section { title: qsTr("Model") - VariablesForm + FactorsForm { - preferredHeight: jaspTheme.smallDefaultVariablesFormHeight - AvailableVariablesList - { - name: "availableTerms" - title: qsTr("Components") - width: parent.width / 4 - source: ['covariates', 'factors'] - } - ModelTermsList { width: parent.width * 5 / 9 } + id: factors + name: "modelTerms" + allowAll: true + nested: nested.checked + allowInteraction: true + initNumberFactors: 2 + baseName: "model" + baseTitle: "Model" + availableVariablesList.source: ['covariates', 'factors'] + startIndex: 0 + availableVariablesListName: "availableTerms" + } + + CheckBox + { + id: nested + label: "Nested" + name: "nested" + checked: true } CheckBox { name: "interceptTerm"; label: qsTr("Include intercept"); checked: true } @@ -74,84 +88,68 @@ Form { title: qsTr("Statistics") + columns: 2 + Group + { + title: qsTr("Model Summary") + CheckBox { name: "rSquaredChange"; label: qsTr("R squared change") } + CheckBox { name: "fChange"; label: qsTr("F change") } + CheckBox { name: "modelAICBIC"; label: qsTr("AIC and BIC") } + CheckBox { name: "residualDurbinWatson"; label: qsTr("Durbin-Watson") } + + } + Group { title: qsTr("Coefficients") - columns: 2 - Layout.columnSpan: 2 - Group + CheckBox { + name: "coefficientEstimate" + label: qsTr("Estimates") + checked: true + onClicked: { if (!checked && bootstrapping.checked) bootstrapping.click() } CheckBox { - name: "coefficientEstimate" - label: qsTr("Estimates") - checked: true - onClicked: { if (!checked && bootstrapping.checked) bootstrapping.click() } - CheckBox + id: bootstrapping + name: "coefficientBootstrap" + label: qsTr("From") + childrenOnSameRow: true + IntegerField { - id: bootstrapping - name: "coefficientBootstrap" - label: qsTr("From") - childrenOnSameRow: true - IntegerField - { - name: "coefficientBootstrapSamples" - defaultValue: 5000 - fieldWidth: 50 - min: 100 - afterLabel: qsTr("bootstraps") - } + name: "coefficientBootstrapSamples" + defaultValue: 5000 + fieldWidth: 50 + min: 100 + afterLabel: qsTr("bootstraps") } } - - CheckBox - { - name: "coefficientCi"; label: qsTr("Confidence intervals") - childrenOnSameRow: true - CIField { name: "coefficientCiLevel" } - } - CheckBox { name: "covarianceMatrix"; label: qsTr("Covariance matrix") } - CheckBox { name: "vovkSellke"; label: qsTr("Vovk-Sellke maximum p-ratio") } } - Group + CheckBox { - CheckBox { name: "modelFit"; label: qsTr("Model fit"); checked: true } - CheckBox { name: "rSquaredChange"; label: qsTr("R squared change") } - CheckBox { name: "descriptives"; label: qsTr("Descriptives") } - CheckBox { name: "partAndPartialCorrelation"; label: qsTr("Part and partial correlations") } - CheckBox { name: "collinearityDiagnostic"; label: qsTr("Collinearity diagnostics") } + name: "coefficientCi"; label: qsTr("Confidence intervals") + childrenOnSameRow: true + CIField { name: "coefficientCiLevel" } } + CheckBox { name: "collinearityStatistic"; label: qsTr("Tolerance and VIF") } + CheckBox { name: "vovkSellke"; label: qsTr("Vovk-Sellke maximum p-ratio") } } Group { - title: qsTr("Residuals") - CheckBox { name: "residualStatistic"; label: qsTr("Statistics") } - CheckBox { name: "residualDurbinWatson"; label: qsTr("Durbin-Watson") } - CheckBox - { - name: "residualCasewiseDiagnostic"; label: qsTr("Casewise diagnostics") - RadioButtonGroup - { - name: "residualCasewiseDiagnosticType" - RadioButton - { - value: "outliersOutside"; label: qsTr("Standard residual >"); checked: true - childrenOnSameRow: true - DoubleField { name: "residualCasewiseDiagnosticZThreshold"; defaultValue: 3 } - } - RadioButton - { - value: "cooksDistance"; label: qsTr("Cook's distance >") - childrenOnSameRow: true - DoubleField { name: "residualCasewiseDiagnosticCooksDistanceThreshold"; defaultValue: 1 } - } - RadioButton { value: "allCases"; label: qsTr("All") } - } - } + title: qsTr("Display") + CheckBox { name: "modelFit"; label: qsTr("Model fit"); checked: true } + CheckBox { name: "descriptives"; label: qsTr("Descriptives") } + CheckBox { name: "partAndPartialCorrelation"; label: qsTr("Part and partial correlations") } + CheckBox { name: "covarianceMatrix"; label: qsTr("Coefficients covariance matrix") } + CheckBox { name: "collinearityDiagnostic"; label: qsTr("Collinearity diagnostics") } + } + + + Common.OutlierComponent { id: outlierComponentt} + } Section diff --git a/inst/qml/RegressionLogistic.qml b/inst/qml/RegressionLogistic.qml index dabde8c8..22be9a43 100644 --- a/inst/qml/RegressionLogistic.qml +++ b/inst/qml/RegressionLogistic.qml @@ -123,7 +123,7 @@ Form CheckBox { name: "vovkSellke"; label: qsTr("Vovk-Sellke maximum p-ratio") } } - CheckBox { name: "multicollinearity"; label: qsTr("Multicollinearity Diagnostics") } + CheckBox { name: "multicollinearity"; label: qsTr("Multicollinearity diagnostics") } } Group diff --git a/inst/qml/common/OutlierComponent.qml b/inst/qml/common/OutlierComponent.qml new file mode 100644 index 00000000..2ec819f6 --- /dev/null +++ b/inst/qml/common/OutlierComponent.qml @@ -0,0 +1,84 @@ +// +// Copyright (C) 2013-2018 University of Amsterdam +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU Affero General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Affero General Public License for more details. +// +// You should have received a copy of the GNU Affero General Public +// License along with this program. If not, see +// . +// + +import QtQuick +import QtQuick.Layouts +import JASP +import JASP.Controls + +// All Analysis forms must be built with the From QML item +Group +{ + title: qsTr("Residuals") + CheckBox + { + name: "residualStatistic" + label: qsTr("Statistics") + visible: analysis === Type.Analysis.LinearRegression + } + + CheckBox + { + name: "residualCasewiseDiagnostic"; label: qsTr("Casewise diagnostics") + columns: 2 + RadioButtonGroup + { + name: "residualCasewiseDiagnosticType" + RadioButton + { + value: "outliersOutside"; label: qsTr("Std. residual >"); checked: true + childrenOnSameRow: true + DoubleField { name: "residualCasewiseDiagnosticZThreshold"; defaultValue: 3 } + } + RadioButton + { + value: "cooksDistance"; label: qsTr("Cook's dist. >") + childrenOnSameRow: true + DoubleField { name: "residualCasewiseDiagnosticCooksDistanceThreshold"; defaultValue: 1 } + } + RadioButton { value: "allCases"; label: qsTr("All") } + } + + Group + { + CheckBox { name: "dfbetas"; label: qsTr("DFBETAS") } + CheckBox { name: "dffits"; label: qsTr("DFFITS") } + CheckBox { name: "covarianceRatio"; label: qsTr("Cov ratio") } + CheckBox { name: "leverage"; label: qsTr("Leverage") } + CheckBox { name: "mahalanobis"; label: qsTr("Mahalanobis") } + + } + + } + + CheckBox + { + id: residualsSavedToData + name: "residualsSavedToData" + text: qsTr("Append residuals to data") + + ComputedColumnField + { + name: "residualsSavedToDataColumn" + text: qsTr("Column name") + placeholderText: qsTr("e.g., residuals") + fieldWidth: 120 + enabled: residualsSavedToData.checked + } + } +} \ No newline at end of file diff --git a/inst/qml/common/Type.qml b/inst/qml/common/Type.qml new file mode 100644 index 00000000..8f0314ef --- /dev/null +++ b/inst/qml/common/Type.qml @@ -0,0 +1,34 @@ +// +// Copyright (C) 2013-2022 University of Amsterdam +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU Affero General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Affero General Public License for more details. +// +// You should have received a copy of the GNU Affero General Public +// License along with this program. If not, see +// . +// + +import QtQuick + +Item +{ + enum Framework + { + Classical, + Bayesian + } + + enum Analysis + { + LinearRegression, + GLM + } +} diff --git a/tests/testthat/_snaps/generalizedlinearmodel/devresqqplot.svg b/tests/testthat/_snaps/generalizedlinearmodel/devresqqplot.svg index 506791e6..b4dc99b9 100644 --- a/tests/testthat/_snaps/generalizedlinearmodel/devresqqplot.svg +++ b/tests/testthat/_snaps/generalizedlinearmodel/devresqqplot.svg @@ -27,55 +27,49 @@ - - - - - - - - - - - - + + + + + + + + + + + + - + --2 --1 -0 -1 -2 +-2 +-1 +0 +1 +2 3 - - - - - + + + + + - - - - - - - + + + + -2 --1.5 --1 --0.5 -0 -0.5 -1 -1.5 -2 +-1 +0 +1 +2 +3 Theoretical quantiles Observed quantiles devResQqPlot diff --git a/tests/testthat/_snaps/generalizedlinearmodel/prsresqqplot.svg b/tests/testthat/_snaps/generalizedlinearmodel/prsresqqplot.svg index 122be629..99a8e4cb 100644 --- a/tests/testthat/_snaps/generalizedlinearmodel/prsresqqplot.svg +++ b/tests/testthat/_snaps/generalizedlinearmodel/prsresqqplot.svg @@ -27,55 +27,49 @@ - - - - - - - - - - - - + + + + + + + + + + + + - + --2 --1 -0 -1 -2 +-2 +-1 +0 +1 +2 3 - - - - - + + + + + - - - - - - - + + + + -2 --1.5 --1 --0.5 -0 -0.5 -1 -1.5 -2 +-1 +0 +1 +2 +3 Theoretical quantiles Observed quantiles prsResQqPlot diff --git a/tests/testthat/_snaps/generalizedlinearmodel/quanresqqplot.svg b/tests/testthat/_snaps/generalizedlinearmodel/quanresqqplot.svg index b73c7c98..fe70eec8 100644 --- a/tests/testthat/_snaps/generalizedlinearmodel/quanresqqplot.svg +++ b/tests/testthat/_snaps/generalizedlinearmodel/quanresqqplot.svg @@ -27,55 +27,49 @@ - - - - - - - - - - - - + + + + + + + + + + + + - + --2 --1 -0 -1 -2 +-2 +-1 +0 +1 +2 3 - - - - - + + + + + - - - - - - - + + + + -2 --1.5 --1 --0.5 -0 -0.5 -1 -1.5 -2 +-1 +0 +1 +2 +3 Theoretical quantiles Observed quantiles quanResQqPlot diff --git a/tests/testthat/_snaps/regressionlinear/field-qq.svg b/tests/testthat/_snaps/regressionlinear/field-qq.svg index 1d87599d..ee0c9dfe 100644 --- a/tests/testthat/_snaps/regressionlinear/field-qq.svg +++ b/tests/testthat/_snaps/regressionlinear/field-qq.svg @@ -27,244 +27,248 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + -3 --2 --1 -0 -1 -2 -3 +-2 +-1 +0 +1 +2 +3 +4 - - - - - - + + + + + + + - - - - - + + + + + + -3 --2 --1 -0 -1 -2 -3 -Theoretical Quantiles -Standardized Residuals +-2 +-1 +0 +1 +2 +3 +4 +Theoretical quantiles +Observed quantiles field-qq diff --git a/tests/testthat/_snaps/regressionlinear/residuals-q-q.svg b/tests/testthat/_snaps/regressionlinear/residuals-q-q.svg index 3481363f..85236306 100644 --- a/tests/testthat/_snaps/regressionlinear/residuals-q-q.svg +++ b/tests/testthat/_snaps/regressionlinear/residuals-q-q.svg @@ -27,144 +27,148 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + -3 --2 --1 -0 -1 -2 -3 +-2 +-1 +0 +1 +2 +3 +4 - - - - - - + + + + + + + - - - - - + + + + + + -3 --2 --1 -0 -1 -2 -3 -Theoretical Quantiles -Standardized Residuals +-2 +-1 +0 +1 +2 +3 +4 +Theoretical quantiles +Observed quantiles residuals-q-q diff --git a/tests/testthat/_snaps/regressionlinear/socialanxiety-qq-plot.svg b/tests/testthat/_snaps/regressionlinear/socialanxiety-qq-plot.svg index 8129066e..f9ca97ae 100644 --- a/tests/testthat/_snaps/regressionlinear/socialanxiety-qq-plot.svg +++ b/tests/testthat/_snaps/regressionlinear/socialanxiety-qq-plot.svg @@ -28,138 +28,138 @@ - - - - - - - - + + + + + + + + - - - - - - - + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - + + + + + + + - + - - - - - + + + + + - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -195,8 +195,8 @@ 1 2 3 -Theoretical Quantiles -Standardized Residuals +Theoretical quantiles +Observed quantiles socialAnxiety-QQ-plot diff --git a/tests/testthat/_snaps/regressionlinearbayesian/qqplot.svg b/tests/testthat/_snaps/regressionlinearbayesian/qqplot.svg index 0a8713c7..e6818583 100644 --- a/tests/testthat/_snaps/regressionlinearbayesian/qqplot.svg +++ b/tests/testthat/_snaps/regressionlinearbayesian/qqplot.svg @@ -27,7 +27,7 @@ - + diff --git a/tests/testthat/test-correlation.R b/tests/testthat/test-correlation.R index fdc09a65..57d5bf30 100644 --- a/tests/testthat/test-correlation.R +++ b/tests/testthat/test-correlation.R @@ -15,6 +15,7 @@ options$scatterPlotStatistic <- TRUE options$sampleSize <- TRUE options$spearman <- TRUE options$effectSize <- TRUE +options$covariance <- TRUE options$variables <- list("contNormal", "contGamma", "contcor1", "debMiss30") set.seed(1) results <- jaspTools::runAnalysis("Correlation", "debug.csv", options) @@ -65,63 +66,66 @@ test_that("Spearman's rho heatmap matches", { test_that("Correlation Table results match", { table <- results[["results"]][["mainTable"]][["data"]] jaspTools::expect_equal_tables(table, - list(-0.0266729903526464, -0.0266666666666667, -0.15041383947394, 0.694237192757787, - 0.0677819401868667, 0.097080506140607, 1, -0.0592696913271387, - -0.0592003859505642, -0.252680329590477, 0.558497687623534, - 0.101534616513362, 0.138832075039338, 1, 100, -0.0341927371158639, - -0.0341794179417942, -0.229059752837501, 0.73526094223706, 0.101287863086583, - 0.163335243866025, 1, "", "", "", + list(-0.0960185736017089, -0.0266729903526464, -0.0266666666666667, + -0.15041383947394, 0.694237192757787, 0.0677819401868667, 0.097080506140607, + 1, -0.0592696913271387, -0.0592003859505643, -0.252680329590477, + 0.558497687623534, 0.101534616513362, 0.138832075039338, 1, + 100, -0.0341927371158639, -0.0341794179417942, -0.229059752837501, + 0.73526094223706, 0.101287863086583, 0.163335243866025, 1, "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - "", "", "", "", 0.0960518800942078, + "", "", "", "", "", + "", "", 0.17245614334619, 0.0960518800942078, 0.0957575757575758, -0.0373731390706183, 0.15805971278439, 0.0671525377474683, 0.22888829058577, 1.26165085338952, 0.16244591532518, 0.161031927910319, - -0.0365419981231702, 0.109479317429059, 0.101534616513362, 0.346490687832583, + -0.03654199812317, 0.109479317429059, 0.101534616513362, 0.346490687832583, 1.51909334263147, 100, 0.143821784353644, 0.142838283828383, -0.0551264633902869, 0.156055917429528, 0.102156998059743, 0.329997969616898, - 1.26907384634445, -0.142995727486937, -0.142028985507246, -0.302753498566225, - 0.0820536231540238, 0.0798902375559328, 0.0186955275517326, - 1.7930869050848, -0.165268370300722, -0.163779936728643, -0.383976976749411, - 0.175488795918533, 0.122169444356305, 0.0740435803355283, 1.20465290217953, - 70, -0.208206182304557, -0.20524888461202, -0.419968595404043, - 0.0883143492445961, 0.1232177656224, 0.0312313683562874, 1.71644871351761, - "contNormal", "", "", "", "", + 1.26907384634445, -4.54234111641073, -0.142995727486937, -0.142028985507246, + -0.302753498566225, 0.0820536231540238, 0.0798902375559328, + 0.0186955275517326, 1.7930869050848, -0.165268370300722, -0.163779936728643, + -0.383976976749411, 0.175488795918533, 0.122169444356305, 0.0740435803355283, + 1.20465290217953, 70, -0.208206182304557, -0.20524888461202, + -0.419968595404043, 0.0883143492445961, 0.1232177656224, 0.0312313683562874, + 1.71644871351761, "contNormal", "", "", "", + "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - "", -0.12919895977801, -0.128484848484848, -0.265695191109086, - 0.0582140897855729, 0.0666069384089473, 0.00872549413938944, - 2.22231002833737, -0.157853482232319, -0.156555303374674, -0.342443190888167, - 0.119832226549265, 0.101534616513362, 0.041127497099264, 1.44695679291394, - 100, -0.185861370097632, -0.183750375037504, -0.3669254548718, - 0.0673279518522942, 0.102488331218907, 0.0131420647686214, 2.02506621791795, - 0.150610965569096, 0.149482401656315, -0.0220394444690113, 0.0672280148907629, + "", "", "", "", "", "", "", "", -0.242747263345942, -0.12919895977801, + -0.128484848484848, -0.265695191109086, 0.0582140897855729, + 0.0666069384089473, 0.00872549413938947, 2.22231002833737, -0.157853482232319, + -0.156555303374674, -0.342443190888167, 0.119832226549265, 0.101534616513362, + 0.0411274970992641, 1.44695679291394, 100, -0.185861370097632, + -0.183750375037504, -0.3669254548718, 0.0673279518522942, 0.102488331218907, + 0.0131420647686214, 2.02506621791795, 6.82842148086829, 0.150610965569096, + 0.149482401656315, -0.0220394444690113, 0.0672280148907629, 0.0796979487434949, 0.321004247781641, 2.02696064848969, 0.173519134850064, - 0.171798366528544, -0.0658332206699671, 0.155001605969274, 0.122169444356305, + 0.171798366528544, -0.065833220669967, 0.155001605969273, 0.122169444356305, 0.39098888887008, 1.27306010334954, 70, 0.214387923136248, 0.211162627941562, -0.0250545433406204, 0.0793767652827101, 0.123275189177231, 0.425046791840888, 1.82929064467251, "contGamma", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - "", "", "", "", "", "", "", "", "", "", "", "", "", + "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - -0.058451570593472, -0.0583850931677019, -0.226151894979294, - 0.474719152682399, 0.0813754587012183, 0.109381708643891, 1, - -0.0906702415181397, -0.0904225863977578, -0.318626758463425, - 0.456613508199801, 0.122169444356305, 0.147689385556226, 1, - 70, -0.102978167463976, -0.10261569416499, -0.329641395147143, - 0.3970672317383, 0.122236141506121, 0.135628607517475, 1, "contcor1", + "", "", -2.14900575670369, -0.058451570593472, + -0.0583850931677019, -0.226151894979294, 0.474719152682399, + 0.0813754587012183, 0.109381708643891, 1, -0.0906702415181398, + -0.0904225863977578, -0.318626758463425, 0.456613508199801, + 0.122169444356305, 0.147689385556226, 1, 70, -0.102978167463976, + -0.10261569416499, -0.329641395147143, 0.3970672317383, 0.122236141506121, + 0.135628607517475, 1, "contcor1", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", - "", "", "", "", "", "", "", "", "", + "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", diff --git a/tests/testthat/test-generalizedlinearmodel.R b/tests/testthat/test-generalizedlinearmodel.R index 8bb0532d..eafc88c2 100644 --- a/tests/testthat/test-generalizedlinearmodel.R +++ b/tests/testthat/test-generalizedlinearmodel.R @@ -14,7 +14,8 @@ addCommonQmlOptions <- function(options) { options, jaspTools:::readQML(file.path(root, "GlmInputComponent.qml")), jaspTools:::readQML(file.path(root, "GlmResidualAnalysisPlotsComponent.qml")), - jaspTools:::readQML(file.path(root, "EmmComponent.qml")) + jaspTools:::readQML(file.path(root, "EmmComponent.qml")), + jaspTools:::readQML(file.path(root, "OutlierComponent.qml")) ) } @@ -94,10 +95,14 @@ options$standardizedResidualOutlierTableTopN <- 3 options$studentizedResidualOutlierTable <- TRUE options$studentizedResidualOutlierTableTopN <- 3 +options$residualCasewiseDiagnostic <- TRUE +options$residualCasewiseDiagnosticType <- "outliersOutside" +options$residualCasewiseDiagnosticZThreshold <- 1 + + options$dfbetas <- TRUE options$dffits <- TRUE options$covarianceRatio <- TRUE -options$cooksDistance <- TRUE options$leverage <- TRUE options$setSeed <- TRUE @@ -224,9 +229,18 @@ test_that("Outlier table based on studentized deviance residuals matches", { test_that("Influential cases table matches", { table <- results[["results"]][["diagnosticsContainer"]][["collection"]][["diagnosticsContainer_influenceTable"]][["data"]] jaspTools::expect_equal_tables(table, - list(4, -0.1795, 0.146, -0.1918, 1.693, 0.0228, 0.2708, - 9, -0.290, 0.750, 1.439, 0.3537, 0.6321, 0.1902, - 10, 0.1632, -0.2761,-0.3971, 1.687, 0.0982, 0.3111))}) + list(-0.607898874134593, 0.562520030505165, 1, 0.0910041875476176, + 0.811136811632276, 0, -0.608545530898244, 0.122014731899716, + 0.0286397537025298, -1.02948417315995, -1.60670300282295, -0.182578461118574, + 0.0432371780205531, 7, 0.0954834805489629, 1.14396823117325, + 0.214285714285714, -0.421414814442793, 0.141203925314399, 0.283760332058426, + -0.341834673218477, -1.10858029927343, 0, -0.290005764938072, + 0.750048649952999, 9, 0.632130380240399, 0.353736621049185, + 0.647058823529412, 1.43918415935164, 0.190243757445277, 0.468418331001346, + 0.717424210957795, 2.32732272718569, 0.484997819985084, -0.707619585051628, + 11, 0.447283433716663, 1.44554646273235, 0.583333333333333, + -0.886667016095175, 0.36304481969478, 0.662151029321139, -0.352325499848782, + -1.2326217135089))}) # multicollinearity table diff --git a/tests/testthat/test-regressionlinear.R b/tests/testthat/test-regressionlinear.R index c6aa20e3..a09927e2 100644 --- a/tests/testthat/test-regressionlinear.R +++ b/tests/testthat/test-regressionlinear.R @@ -7,16 +7,33 @@ context("Linear Regression") # - stepwise methods (currently gives an error if I set p entry too high) # - plots handle errors -test_that("Main table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") +initOptsLinReg <- function() { + options <- jaspTools::analysisOptions("RegressionLinear") + options$dependent <- "contNormal" options$covariates <- "contGamma" - options$weights <- "facFifty" options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) + list(components= NULL, name="model0", title = "Model 0"), + list(components=list("contGamma"), name="model1", title = "Model 1") ) + options$rSquaredChange <- TRUE + options$fChange <- TRUE + options$residualDurbinWatson <- FALSE + options$residualCasewiseDiagnostic <- FALSE + options$residualsSavedToData <- FALSE + options$residualsSavedToDataColumn <- FALSE + options$residualStatistic <- FALSE + + return(options) +} + +test_that("Main table results match", { + options <- initOptsLinReg() + + options$weights <- "facFifty" options$residualDurbinWatson <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_summaryTable"]][["data"]] jaspTools::expect_equal_tables(table, @@ -29,17 +46,15 @@ test_that("Main table results match", { }) test_that("Coefficients table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$coefficientEstimate <- TRUE options$coefficientCi <- TRUE options$coefficientCiLevel <- 0.9 options$collinearityDiagnostic <- TRUE + options$collinearityStatistic <- TRUE options$vovkSellke <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_coeffTable"]][["data"]] jaspTools::expect_equal_tables(table, @@ -55,14 +70,12 @@ test_that("Coefficients table results match", { }) test_that("ANOVA table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() + options$dependent <- "debCollin1" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) options$modelFit <- TRUE options$vovkSellke <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_anovaTable"]][["data"]] jaspTools::expect_equal_tables(table, @@ -73,32 +86,28 @@ test_that("ANOVA table results match", { }) test_that("Coefficients Covariance table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" + options <- initOptsLinReg() options$covariates <- c("contGamma", "contcor1") options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE), - list(components="contcor1", isNuisance=FALSE) + list(components=list("contGamma"), name="model0", title = "Model 0"), + list(components=list("contGamma", "contcor1"), name="model1", title = "Model 1") ) options$covarianceMatrix <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_coeffCovMatrixTable"]][["data"]] jaspTools::expect_equal_tables(table, - list("TRUE", 0.00490486111017858, 0.00116294327838645, "H", - "contGamma", "FALSE", "", 0.0112500585702943, "H", - "contcor1") + list("FALSE", 0.00485075592634222, "", "M", "contGamma", "TRUE", + 0.00490486111017858, 0.00116294327838645, "M", "contGamma", + "FALSE", "", 0.0112500585702943, "M", "contcor1") ) }) test_that("Descriptive table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$descriptives <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_descriptivesTable"]][["data"]] jaspTools::expect_equal_tables(table, @@ -108,28 +117,32 @@ test_that("Descriptive table results match", { }) test_that("Part and Partial Correlations table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" + options <- initOptsLinReg() + options$covariates <- c("debCollin2", "contGamma") options$modelTerms <- list( - list(components="debCollin2", isNuisance=FALSE), - list(components="contGamma", isNuisance=FALSE) + list(components=list("debCollin2"), name="model0", title = "Model 0"), + list(components=list("contGamma", "debCollin2"), name="model1", title = "Model 1") ) options$partAndPartialCorrelation <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_partialCorTable"]][["data"]] jaspTools::expect_equal_tables(table, - list("TRUE", "H", "debCollin2", -0.0198338559459939, -0.0198687032851428, - "FALSE", "H", "contGamma", -0.0617145529439823, -0.0617173111983368) + list("FALSE", "M", "debCollin2", -0.0094541786161782, -0.0094541786161782, + "TRUE", "M", "contGamma", -0.0617145529439822, -0.0617173111983367, + "FALSE", "M", "debCollin2", -0.0198338559459935, -0.0198687032851424 + ) ) }) test_that("Collinearity Diagonistic table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" + options <- initOptsLinReg() + options$covariates <- "contcor1" options$modelTerms <- list( - list(components="contcor1", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components="contcor1", name="model1", title = "Model 1") ) options$collinearityDiagnostic <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) @@ -142,11 +155,12 @@ test_that("Collinearity Diagonistic table results match", { }) test_that("Residuals Statistics table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "contNormal" options$covariates <- "contcor1" options$modelTerms <- list( - list(components="contcor1", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components="contcor1", name="model1", title = "Model 1") ) options$residualStatistic <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) @@ -162,32 +176,30 @@ test_that("Residuals Statistics table results match", { }) test_that("Casewise Diagnostics table results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "contNormal" options$covariates <- "contOutlier" options$modelTerms <- list( - list(components="contOutlier", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components="contOutlier", name="model1", title = "Model 1") ) options$residualCasewiseDiagnostic <- TRUE options$residualCasewiseDiagnosticType <- "outliersOutside" options$residualCasewiseDiagnosticZThreshold <- 3 results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) - table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_casewiseTable"]][["data"]] + table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_influenceTable"]][["data"]] jaspTools::expect_equal_tables(table, - list(55, 3.34810934796608, 3.356094448, -0.187237305306683, 3.54333175330668, - 0.0577123260439598, 83, 3.22377600371253, 2.958797116, -0.143545494526366, - 3.10234261052637, 1.15289593050314) + list(55, 0.0577123260439598, 3.356094448, -0.187237305306683, 3.54333175330668, + 3.34810934796609, 0, 83, 1.15289593050314, 2.958797116, -0.143545494526367, + 3.10234261052637, 3.22377600371253) ) }) test_that("Residuals vs. Dependent plot matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$residualVsDependentPlot <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) testPlot <- results[["state"]][["figures"]][[1]][["obj"]] @@ -195,52 +207,39 @@ test_that("Residuals vs. Dependent plot matches", { }) test_that("Residuals vs. Covariates plot matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$residualVsCovariatePlot <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) testPlot <- results[["state"]][["figures"]][[1]][["obj"]] jaspTools::expect_equal_plots(testPlot, "residuals-covariates") }) test_that("Residuals vs. Predicted plot matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$residualVsFittedPlot <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) testPlot <- results[["state"]][["figures"]][[1]][["obj"]] jaspTools::expect_equal_plots(testPlot, "residuals-predicted") }) test_that("Standardized Residuals Histogram matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$residualHistogramPlot <- TRUE options$residualHistogramStandardizedPlot <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) testPlot <- results[["state"]][["figures"]][[1]][["obj"]] jaspTools::expect_equal_plots(testPlot, "residuals-histogram") }) test_that("Q-Q Plot Standardized Residuals matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" - options$covariates <- "contGamma" - options$modelTerms <- list( - list(components="contGamma", isNuisance=FALSE) - ) + options <- initOptsLinReg() + options$residualQqPlot <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) testPlot <- results[["state"]][["figures"]][[1]][["obj"]] @@ -248,11 +247,12 @@ test_that("Q-Q Plot Standardized Residuals matches", { }) test_that("Marginal effects plot matches", { - options <- jaspTools::analysisOptions("RegressionLinear") - options$dependent <- "contNormal" + options <- initOptsLinReg() + options$covariates <- "contcor1" options$modelTerms <- list( - list(components="contcor1", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components="contcor1", name="model1", title = "Model 1") ) options$marginalPlot <- TRUE options$marginalPlotCi <- TRUE @@ -266,63 +266,91 @@ test_that("Marginal effects plot matches", { }) test_that("Analysis handles errors", { - options <- jaspTools::analysisOptions("RegressionLinear") - + options <- initOptsLinReg() + options$dependent <- "debInf" options$covariates <- "contGamma" - options$modelTerms <- list(list(components="contGamma", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components="contGamma", name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Inf dependent check") options$dependent <- "contNormal" options$covariates <- "debInf" - options$modelTerms <- list(list(components="debInf", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) + options$modelTerms <- list(list(components="debInf", name="model0", title = "Model 0")) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Inf covariate check") options$covariates <- "contGamma" options$weights <- "debInf" - options$modelTerms <- list(list(components="contGamma", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Inf weights check") options$dependent <- "debSame" options$covariates <- "contGamma" options$weights <- "" - options$modelTerms <- list(list(components="contGamma", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="No variance dependent check") options$dependent <- "contNormal" options$covariates <- "debSame" - options$modelTerms <- list(list(components="debSame", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="No variance covariate check") options$dependent <- "contGamma" options$covariates <- "contcor1" options$weights <- "contNormal" - options$modelTerms <- list(list(components="contcor1", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Negative weights check") options$dependent <- "debNaN" options$covariates <- "contcor1" options$weights <- "" - options$modelTerms <- list(list(components="contcor1", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Too few obs dependent check") options$dependent <- "contGamma" options$covariates <- "debNaN" - options$modelTerms <- list(list(components="debNaN", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Too few obs covariate check") options$dependent <- "contGamma" options$covariates <- "contNormal" options$weights <- "debNaN" - options$modelTerms <- list(list(components="contNormal", isNuisance=FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list(options$covariates), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) expect_identical(results[["status"]], "validationError", label="Too few obs weights check") @@ -330,24 +358,26 @@ test_that("Analysis handles errors", { options$dependent <- "contNormal" options$covariates <- c("debCollin2", "debCollin3") options$modelTerms <- list( - list(components="debCollin2", isNuisance=FALSE), - list(components="debCollin3", isNuisance=FALSE) - ) + list(components=list("debCollin2"), name="model0", title = "Model 0"), + list(components=list("debCollin2", "debCollin3"), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) results$results$errorMessage expect_identical(results[["status"]], "validationError", label="Perfect correlation check") }) test_that("Analysis handles categorical predictors in model summary table", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "contNormal" options$covariates <- "contcor1" options$factors <- "facFive" options$modelTerms <- list( - list(components="contcor1", isNuisance=TRUE), - list(components="facFive", isNuisance=FALSE) - ) + list(components=list("contcor1"), name="model0", title = "Model 0"), + list(components=list("contcor1", "facFive"), name="model1", title = "Model 1") + ) options$rSquaredChange <- TRUE + options$fChange <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", "test.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_summaryTable"]][["data"]] @@ -363,34 +393,34 @@ test_that("Analysis handles categorical predictors in model summary table", { test_that("Part And Partial Correlations table results match", { # Part and partial correlations, including categorical predictors, verified with SPSS, # see pdf doc in https://github.com/jasp-stats/jasp-issues/issues/1638 - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$covariates <- c("education", "prestige") options$dependent <- "income" options$factors <- "occup_type" - options$modelTerms <- list(list(components = "education", isNuisance = FALSE), list( - components = "prestige", isNuisance = FALSE), list(components = "occup_type", - isNuisance = FALSE)) + options$modelTerms <- list( + list(components=list("education", "prestige", "occup_type"), name="model0", title = "Model 0") + ) options$partAndPartialCorrelation <- TRUE set.seed(1) results <- jaspTools::runAnalysis("RegressionLinear", "Duncan.csv", options) table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_partialCorTable"]][["data"]] jaspTools::expect_equal_tables(table, - list("TRUE", "H", "education", -0.107546637580272, -0.218295377900754, - "FALSE", "H", "prestige", 0.363479641784283, 0.603066145085564, - "FALSE", "H", "occup_type", 0.236073658532065, 0.440752114010841 + list("FALSE", "M", "education", -0.107546637580271, -0.218295377900751, + "FALSE", "M", "prestige", 0.363479641784282, 0.603066145085563, + "FALSE", "M", "occup_type", 0.236073658532064, 0.44075211401084 )) }) test_that("Bootstrapping runs", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "contNormal" options$covariates <- "contcor1" options$factors <- "facFive" options$modelTerms <- list( - list(components="contcor1", isNuisance=TRUE), - list(components="facFive", isNuisance=FALSE) - ) + list(components=list("contcor1"), name="model0", title = "Model 0"), + list(components=list("contcor1", "facFive"), name="model1", title = "Model 1") + ) options$coefficientBootstrap <- TRUE options$coefficientBootstrapSamples <- 100 options$coefficientCi <- TRUE @@ -423,16 +453,17 @@ test_that("Bootstrapping runs", { }) test_that("Marginal effects plots works with interactions", { - options <- analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$.meta <- list(covariates = list(shouldEncode = TRUE), dependent = list( shouldEncode = TRUE), factors = list(shouldEncode = TRUE), modelTerms = list(shouldEncode = TRUE), weights = list(shouldEncode = TRUE)) options$covariates <- c("contGamma", "contcor1") options$dependent <- "contNormal" options$marginalPlot <- TRUE - options$modelTerms <- list(list(components = "contGamma", isNuisance = FALSE), list( - components = "contcor1", isNuisance = FALSE), list(components = c("contGamma", - "contcor1"), isNuisance = FALSE)) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("contGamma", "contcor1", list("contGamma", "contcor1")), name="model1", title = "Model 1") + ) set.seed(1) results <- runAnalysis("RegressionLinear", "test.csv", options) @@ -451,12 +482,16 @@ test_that("Marginal effects plots works with interactions", { # Chapter 1 test_that("Fields Book - Chapter 1 results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "sales" options$covariates <- "adverts" options$modelTerms <- list( - list(components="adverts", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("adverts"), name="model1", title = "Model 1") ) + options$fChange <- FALSE + options$rSquaredChange <- FALSE + results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Album Sales.csv", options) output1 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_summaryTable"]][["data"]] jaspTools::expect_equal_tables(output1, @@ -482,10 +517,9 @@ test_that("Fields Book - Chapter 1 results match", { options$covariates <- c("adverts", "airplay", "attract") options$modelTerms <- list( - list(components="adverts", isNuisance=TRUE), - list(components="airplay", isNuisance=FALSE), - list(components="attract", isNuisance=FALSE) - ) + list(components=list("adverts"), name="model0", title = "Model 0"), + list(components=list("adverts", "airplay", "attract"), name="model1", title = "Model 1") + ) results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Album Sales.csv", options) output4 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_summaryTable"]][["data"]] jaspTools::expect_equal_tables(output4, @@ -521,15 +555,16 @@ test_that("Fields Book - Chapter 1 results match", { # Chapter 2 test_that("Fields Book - Chapter 2 results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "sales" options$covariates <- c("adverts", "airplay", "attract") options$modelTerms <- list( - list(components="adverts", isNuisance=TRUE), - list(components="airplay", isNuisance=FALSE), - list(components="attract", isNuisance=FALSE) - ) + list(components=list("adverts"), name="model0", title = "Model 0"), + list(components=list("adverts", "airplay", "attract"), name="model1", title = "Model 1") + ) options$rSquaredChange <- TRUE + options$fChange <- TRUE + options$coefficientCi <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Album Sales.csv", options) output4 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_summaryTable"]][["data"]] @@ -551,14 +586,13 @@ test_that("Fields Book - Chapter 2 results match", { # Chapter 3 test_that("Fields Book - Chapter 3 results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "sales" options$covariates <- c("adverts", "airplay", "attract") options$modelTerms <- list( - list(components="adverts", isNuisance=TRUE), - list(components="airplay", isNuisance=FALSE), - list(components="attract", isNuisance=FALSE) - ) + list(components=list("adverts"), name="model0", title = "Model 0"), + list(components=list("adverts", "airplay", "attract"), name="model1", title = "Model 1") + ) options$residualVsFittedPlot <- TRUE options$partialResidualPlot <- TRUE options$residualHistogramPlot <- TRUE @@ -568,6 +602,7 @@ test_that("Fields Book - Chapter 3 results match", { options$residualCasewiseDiagnosticType <- "outliersOutside" options$residualCasewiseDiagnosticZThreshold <- 2 options$coefficientCi <- TRUE + results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Album Sales.csv", options) figure3 <- results[["state"]][["figures"]][[1]][["obj"]] # Residuals vs. Predicted jaspTools::expect_equal_plots(figure3, "field-residuals-predicted") @@ -581,7 +616,7 @@ test_that("Fields Book - Chapter 3 results match", { jaspTools::expect_equal_plots(figure5a, "field-residuals-histogram") figure5b <- results[["state"]][["figures"]][[3]][["obj"]] # Q-Q-Plot jaspTools::expect_equal_plots(figure5b, "field-qq") - output1 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_casewiseTable"]][["data"]] + output1 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_influenceTable"]][["data"]] jaspTools::expect_equal_tables(output1, list(1, 2.177404, 330, 229.9203, 100.0797, 0.05870388, 2, -2.323083, 120, 228.949, -108.949, 0.01088943, @@ -597,14 +632,13 @@ test_that("Fields Book - Chapter 3 results match", { 200, -2.088044, 110, 207.2061, -97.20606, 0.02513455) ) - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "sales" options$covariates <- c("adverts", "airplay", "attract") options$modelTerms <- list( - list(components="adverts", isNuisance=TRUE), - list(components="airplay", isNuisance=FALSE), - list(components="attract", isNuisance=FALSE) - ) + list(components=list("adverts"), name="model0", title = "Model 0"), + list(components=list("adverts", "airplay", "attract"), name="model1", title = "Model 1") + ) options$residualCasewiseDiagnostic <- TRUE options$residualCasewiseDiagnosticType <- "allCases" options$coefficientBootstrap <- TRUE @@ -612,7 +646,7 @@ test_that("Fields Book - Chapter 3 results match", { set.seed(1) # For Bootstrapping Unit Tests options$coefficientCi <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Album Sales.csv", options) - figure10 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_casewiseTable"]][["data"]] + figure10 <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_influenceTable"]][["data"]] figure10 <- list(figure10[[1]]$cooksD, figure10[[2]]$cooksD, figure10[[3]]$cooksD, figure10[[4]]$cooksD, figure10[[5]]$cooksD, figure10[[6]]$cooksD, figure10[[7]]$cooksD, figure10[[8]]$cooksD, figure10[[9]]$cooksD, figure10[[10]]$cooksD, figure10[[11]]$cooksD, figure10[[12]]$cooksD, @@ -637,13 +671,17 @@ test_that("Fields Book - Chapter 3 results match", { # "", "attract", -0.2110928, 11.08634, 2.234141, 6.502079, 15.13605) # ) - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "spai" options$covariates <- c("tosca", "obq") options$modelTerms <- list( - list(components="tosca", isNuisance=TRUE), - list(components="obq", isNuisance=FALSE) + list(components="tosca", name="model0", title = "Model 0"), + list(components="obq", name="model0", title = "Model 0") ) + options$modelTerms <- list( + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("tosca", "obq"), name="model1", title = "Model 1") + ) options$residualVsFittedPlot <- TRUE options$partialResidualPlot <- TRUE options$residualQqPlot <- TRUE @@ -662,13 +700,13 @@ test_that("Fields Book - Chapter 3 results match", { # Chapter 4 test_that("Fields Book - Chapter 4 results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "Happiness" options$covariates <- c("dummy1", "dummy2") options$modelTerms <- list( - list(components="dummy1", isNuisance=FALSE), - list(components="dummy2", isNuisance=FALSE) - ) + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("dummy1", "dummy2"), name="model1", title = "Model 1") + ) options$coefficientCi <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Puppies Dummy.csv", options) output1a <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_anovaTable"]][["data"]] @@ -690,13 +728,13 @@ test_that("Fields Book - Chapter 4 results match", { # Chapter 5 test_that("Fields Book - Chapter 5 results match", { - options <- jaspTools::analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$dependent <- "Happiness" options$covariates <- c("Dummy1", "Dummy2") options$modelTerms <- list( - list(components="Dummy1", isNuisance=FALSE), - list(components="Dummy2", isNuisance=FALSE) - ) + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("Dummy1", "Dummy2"), name="model1", title = "Model 1") + ) options$coefficientCi <- TRUE results <- jaspTools::runAnalysis("RegressionLinear", dataset = "Puppies Contrast.csv", options) output1a <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_anovaTable"]][["data"]] @@ -719,20 +757,19 @@ test_that("VIF is correct when the model contains factors", { # test issue reported in https://forum.cogsci.nl/discussion/comment/27675 # previously, we computed the variance inflation factor (vif) in an incorrect manner when categorical variables were in play. - options <- analysisOptions("RegressionLinear") + options <- initOptsLinReg() options$collinearityDiagnostic <- TRUE + options$collinearityStatistic <- TRUE options$covariates <- c("contcor1", "contcor2", "contGamma") options$dependent <- "contNormal" options$factors <- c("contBinom", "facFive") options$modelTerms <- list( - list(components = "contcor1", isNuisance = FALSE), - list(components = "contcor2", isNuisance = FALSE), - list(components = "contGamma", isNuisance = FALSE), - list(components = "contBinom", isNuisance = FALSE), - list(components = "facFive", isNuisance = FALSE), - list(components = c("contBinom", "facFive"), isNuisance = FALSE), - list(components = c("contcor1", "facFive"), isNuisance = FALSE) - ) + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("contcor1", "contcor2", "contGamma", "contBinom", "facFive", + list("contBinom", "facFive"), + list("contcor1", "facFive")), + name="model1", title = "Model 1") + ) set.seed(1) results <- runAnalysis("RegressionLinear", "debug.csv", options) diff --git a/tests/testthat/test-regressionlinearVerification.R b/tests/testthat/test-regressionlinearVerification.R index 0a7298eb..7d91c8f5 100644 --- a/tests/testthat/test-regressionlinearVerification.R +++ b/tests/testthat/test-regressionlinearVerification.R @@ -14,9 +14,15 @@ options$dependent <- "Y" options$covariates <- "MeanCenteredX" options$modelTerms <- list( - list(components="MeanCenteredX", isNuisance=FALSE) + list(components=NULL, name="model0", title = "Model 0"), + list(components=list("MeanCenteredX"), name="model1", title = "Model 1") ) options$descriptives <- TRUE +options$residualDurbinWatson <- FALSE +options$residualCasewiseDiagnostic <- FALSE +options$residualsSavedToData <- FALSE +options$residualsSavedToDataColumn <- FALSE +options$residualStatistic <- FALSE results <- jaspTools::runAnalysis("RegressionLinear", "Regression.csv", options)