Allow constant variance

Author: kvelleby

R/formula.R

@@ -12,22 +12,22 @@ parse_hetero_formula <- function(formula) {
 	f_char <- Reduce(paste, deparse(formula))
 
 
-	# Check for | separator
-
+	# Check for | separator - if missing, default to intercept-only variance
 	if (!grepl("\\|", f_char)) {
-		stop("Formula must contain '|' separating mean and variance equations.\n",
-				 "  Usage: y ~ x1 + x2 | z1 + z2", call. = FALSE)
-	}
+		# No variance equation specified - use intercept-only (homoscedastic)
+		mean_part <- f_char
+		var_part <- "1"
+	} else {
+		# Split: everything before | is mean equation, after is variance equation
+		parts <- strsplit(f_char, "\\|")[[1]]
+		if (length(parts) != 2) {
+			stop("Formula must have exactly one '|' separator.", call. = FALSE)
+		}
 
-	# Split: everything before | is mean equation, after is variance equation
-	parts <- strsplit(f_char, "\\|")[[1]]
-	if (length(parts) != 2) {
-		stop("Formula must have exactly one '|' separator.", call. = FALSE)
+		mean_part <- trimws(parts[1])  # "y ~ x1 + x2"
+		var_part <- trimws(parts[2])   # "z1 + z2"
 	}
 
-	mean_part <- trimws(parts[1])  # "y ~ x1 + x2"
-	var_part <- trimws(parts[2])   # "z1 + z2"
-
 	# Parse mean equation (standard formula)
 	mean_formula <- as.formula(mean_part)
 	lhs <- all.vars(mean_formula[[2]])  # Response variable

R/methods.R

@@ -109,9 +109,17 @@ summary.hetero_fit <- function(object, ...) {
 	}
 
 	# Heteroscedasticity assessment: SD of Pearson residuals by quantile of fitted variance
-	var_quantiles <- cut(pred$sigma, breaks = quantile(pred$sigma, probs = c(0, 0.25, 0.5, 0.75, 1)),
-											 include.lowest = TRUE, labels = c("Q1", "Q2", "Q3", "Q4"))
-	resid_sd_by_var <- tapply(pearson_resid, var_quantiles, sd)
+	# Handle constant sigma case (homoscedastic model with intercept-only variance)
+	sigma_range <- diff(range(pred$sigma))
+	if (sigma_range < .Machine$double.eps * 100) {
+		# Constant sigma - no heteroscedasticity to assess
+		resid_sd_by_var <- c(Q1 = sd(pearson_resid))
+		names(resid_sd_by_var) <- "All (constant variance)"
+	} else {
+		var_quantiles <- cut(pred$sigma, breaks = quantile(pred$sigma, probs = c(0, 0.25, 0.5, 0.75, 1)),
+												 include.lowest = TRUE, labels = c("Q1", "Q2", "Q3", "Q4"))
+		resid_sd_by_var <- tapply(pearson_resid, var_quantiles, sd)
+	}
 
 	# Build summary object
 	result <- list(