Skip to content

Commit

Permalink
Release fixes (#277)
Browse files Browse the repository at this point in the history 
* a long model formula does not break checking whether a model is intercept-only

* empty model does not crash .vif.default but instead gives an informative error msg

* failing wald test does not crash the entire coefficients table

* failing vif does not crash analysis

* failing durbin-watson does not crash analysis

* fix typo
  • Loading branch information
@Kucharssim
Kucharssim authored Dec 12, 2023
1 parent bb5d60b commit 50e718f
Show file tree
Hide file tree
Showing 4 changed files with 59 additions and 46 deletions.
51 changes: 18 additions & 33 deletions R/commonglm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,8 +155,7 @@

.mcFadden <- function(glmModel, nullModel) {
# https://eml.berkeley.edu/reprints/mcfadden/zarembka.pdf
rightSide <- deparse(glmModel[["formula"]][[3]])
if (length(rightSide == 1) && rightSide %in% c("1", "0")) {
if (.isInterceptOnly(glmModel)) {
# intercept-only model needs fix because of computer precision limits
return(0)
} else
Expand All @@ -165,8 +164,7 @@

.nagelkerke <- function(glmModel, nullModel) {
# https://doi.org/10.1093/biomet/78.3.691
rightSide <- deparse(glmModel[["formula"]][[3]])
if (length(rightSide == 1) && rightSide %in% c("1", "0")) {
if (.isInterceptOnly(glmModel)) {
# intercept-only model needs fix because of computer precision limits
return(NULL)
} else {
Expand All @@ -187,8 +185,7 @@
}

.coxSnell <- function(glmModel, nullModel) {
rightSide <- deparse(glmModel[["formula"]][[3]])
if (length(rightSide == 1) && rightSide %in% c("1", "0")) {
if (.isInterceptOnly(glmModel)) {
# intercept-only model needs fix because of computer precision limits
return(NULL)
} else {
Expand Down Expand Up @@ -559,32 +556,20 @@
dw
}

# multicollineary statistics, also from package car 3.0-12
.vif.default <- function(mod, ...) {
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)
# test if a model is intercept only
.isInterceptOnly <- function(model) {
terms <- try(terms(model))
if (jaspBase::isTryError(terms)) {
terms <- try(terms(model[["formula"]]))

if (jaspBase::isTryError(terms)) {
stop("Something went wrong; Cannot tell if a model is intercept-only.")
}
}
if (all(result[, 2] == 1)) result <- result[, 1]
else result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
result

intercept <- attr(terms, "intercept")
labels <- attr(terms, "term.labels")

# test if we have intercept and no other predictors
return(intercept && length(labels) == 0)
}
5 changes: 5 additions & 0 deletions R/glmCommonFunctions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,6 +229,11 @@
# 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)
Expand Down
20 changes: 16 additions & 4 deletions R/regressionlinear.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1077,7 +1077,12 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) {
durbinWatson <- list(r = NaN, dw = NaN, p = NaN)

if (!is.null(fit)) {
durbinWatson <- .durbinWatsonTest.lm(fit, alternative = c("two.sided"))
# TODO: Make some nicer error messsage/footnote when durbin watson computation fails
durbinWatson <- try(.durbinWatsonTest.lm(fit, alternative = c("two.sided")))

if (jaspBase::isTryError(durbinWatson)) {
return(list(r = NaN, dw = NaN, p = NaN))
}

if (weights == "") # if regression is not weighted, calculate p-value with lmtest (car method is unstable)
# TODO: this can fail when there are many interactions between factors. Do we want to show a footnote about that?
Expand Down Expand Up @@ -1233,9 +1238,9 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) {

}

# get translation for (Intercept)
# get translation for (Intercept)
name <- names[i]
if (identical(name, "(Intercept)"))
if (identical(name, "(Intercept)"))
name <- gettext("(Intercept)")

row <- list(
Expand Down Expand Up @@ -1333,7 +1338,14 @@ RegressionLinearInternal <- function(jaspResults, dataset = NULL, options) {
}

# we can actually compute things
result <- car::vif(fit)
result <- try(car::vif(fit))

if (jaspBase::isTryError(result)) {
nas <- rep(NA, noTerms)
names(nas) <- labels(terms(fit))
result <- list(VIF = nas, tolerance = nas)
return(result)
}

VIF <- if (is.matrix(result)) {
result[, 3L]
Expand Down
29 changes: 20 additions & 9 deletions R/regressionlogistic.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,6 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
modelSummary$addColumnInfo(name = "cas", title = gettextf("Cox & Snell R%s","\u00B2"), type = "number")

jaspResults[["modelSummary"]] <- modelSummary

res <- try(.reglogisticModelSummaryFill(jaspResults, dataset, options, ready))

.reglogisticSetError(res, modelSummary)
Expand Down Expand Up @@ -506,7 +505,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
s[3] <- s[1]/s[2] # new z
s[4] <- 2*pnorm(-abs(s[3])) # new p
}
waldtest <- mdscore::wald.test(glmObj[[2]], term = 1)
waldtest <- .waldTest(glmObj[[2]], term = 1)
jaspResults[["estimatesTable"]]$addRows(
list( param = .formatTerm(rn, glmObj[[2]]),
est = s[1],
Expand All @@ -515,11 +514,12 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
or = exp(s[1]),
zval = s[3],
pval = s[4],
waldsta = as.numeric(waldtest$W),
waldsta = as.numeric(waldtest),
walddf = as.numeric(1),
vsmpr = VovkSellkeMPR(s[4]),
cilo = expon(s[1] - alpha * s[2]),
ciup = expon(s[1] + alpha * s[2])))

} else {
if (options$robustSe) {
s[,2] <- unname(.glmRobustSE(glmObj[[2]])) # new se
Expand All @@ -528,7 +528,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
}
for (i in seq_along(rn)) {

waldtest <- mdscore::wald.test(glmObj[[2]], term = i)
waldtest <- .waldTest(glmObj[[2]], term = i)

jaspResults[["estimatesTable"]]$addRows(list(
param = .formatTerm(rn[i], glmObj[[2]]),
Expand All @@ -538,7 +538,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
or = exp(s[i,1]),
zval = s[i,3],
pval = s[i,4],
waldsta = as.numeric(waldtest$W),
waldsta = as.numeric(waldtest),
walddf = as.numeric(1),
vsmpr = VovkSellkeMPR(s[i,4]),
cilo = expon(s[i,1] - alpha * s[i,2]),
Expand Down Expand Up @@ -567,7 +567,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
s[3] <- s[1]/s[2] # new z
s[4] <- 2*pnorm(-abs(s[3])) # new p
}
waldtest <- mdscore::wald.test(mObj, term = 1)
waldtest <- .waldTest(mObj, term = 1)

jaspResults[["estimatesTable"]]$addRows(list(
model = as.character(midx),
Expand All @@ -578,7 +578,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
or = exp(s[1]),
zval = s[3],
pval = s[4],
waldsta = as.numeric(waldtest$W),
waldsta = as.numeric(waldtest),
walddf = as.numeric(1),
vsmpr = VovkSellkeMPR(s[4]),
cilo = expon(s[1] - alpha * s[2]),
Expand All @@ -593,7 +593,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
}
for (i in seq_along(rn)) {

waldtest <- mdscore::wald.test(mObj, term = i)
waldtest <- .waldTest(mObj, term = i)
jaspResults[["estimatesTable"]]$addRows(list(
model = as.character(midx),
param = .formatTerm(rn[i], mObj),
Expand All @@ -603,7 +603,7 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
or = exp(s[i,1]),
zval = s[i,3],
pval = s[i,4],
waldsta = as.numeric(waldtest$W),
waldsta = as.numeric(waldtest),
walddf = as.numeric(1),
vsmpr = VovkSellkeMPR(s[i,4]),
cilo = expon(s[i,1] - alpha * s[i,2]),
Expand Down Expand Up @@ -1382,3 +1382,14 @@ RegressionLogisticInternal <- function(jaspResults, dataset = NULL, options, ...
}
return(list(ciTitle = ciTitle, seTitle = seTitle, multimod = multimod, paramtitle = paramtitle))
}

# wrapper for wald test
.waldTest <- function(...) {
result <- try(mdscore::wald.test(...))
if (jaspBase::isTryError(result)) {
result <- NA
return(result)
}

return(result[["W"]])
}

0 comments on commit 50e718f

Please sign in to comment.