feat(stats): elastic net and lasso via coordinate descent
- ZQElasticNet: cyclic coordinate descent with residual updates, glmnet-convention objective ((1/2n)RSS + lambda(alpha*l1 + (1-alpha)/2*l2)), internal predictor standardization (1/n variance), unpenalized intercept, coefficients reported on the original scale; lambdaMax helper - matched glmnet's gaussian y-standardization quirk deliberately: the L1 penalty is invariant to it but the effective ridge penalty scales by 1/sd(y) — without this, alpha<1 fits diverge from glmnet by ~10% - engine: 'elasticnet y x…, lambda(#) [alpha(#)]' and 'lasso' (alpha fixed at 1); missing lambda() errors with the data's lambda_max as a hint; predict works afterwards, margins refuses (no VCE) - fixtures: glmnet 5.0 at thresh 1e-15 over deliberately correlated regressors; coefficients match at 1e-6 (documented tolerance for penalized iterative solvers) and the selection pattern (which coefficients are exactly zero) matches exactly - ZQCoefficient gains a public initializer - 104 tests green (swift test and xcodebuild with GPU suites) Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Showing 9 changed files with +482 and −0
modified
MetrikaKit/Sources/ZQEngine/Session.swift
+109 −0
@@ -196,6 +196,8 @@ public actor ZQSession { | ||
| 196 | 196 | case "ivregress": return try handleIVRegress(command) |
| 197 | 197 | case "predict": return try handlePredict(command) |
| 198 | 198 | case "margins": return try handleMargins(command) |
| 199 | + case "elasticnet": return try handleElasticNet(command, defaultAlpha: nil) | |
| 200 | + case "lasso": return try handleElasticNet(command, defaultAlpha: 1) | |
| 199 | 201 | case "logit": return try handleGLM(command, family: .logit) |
| 200 | 202 | case "probit": return try handleGLM(command, family: .probit) |
| 201 | 203 | case "poisson": return try handleGLM(command, family: .poisson) |
@@ -1559,6 +1561,113 @@ public actor ZQSession { | ||
| 1559 | 1561 | return ZQResult(text: note, scalars: ["N": Double(n - missingCount)]) |
| 1560 | 1562 | } |
| 1561 | 1563 | |
| 1564 | + /// `elasticnet y x…, lambda(#) [alpha(#)]` and `lasso y x…, lambda(#)` | |
| 1565 | + /// — penalized linear regression, glmnet conventions. Without | |
| 1566 | + /// lambda(), reports λ_max (the smallest all-zero penalty) as a hint. | |
| 1567 | + private func handleElasticNet( | |
| 1568 | + _ command: ZQCommand, defaultAlpha: Double? | |
| 1569 | + ) throws -> ZQResult { | |
| 1570 | + let sample = try buildRegressionSample(command, clusterVariable: nil) | |
| 1571 | + | |
| 1572 | + var alpha = defaultAlpha ?? 1 | |
| 1573 | + if let text = command.option("alpha")?.firstArgument { | |
| 1574 | + guard defaultAlpha == nil else { | |
| 1575 | + throw ZQEngineError("lasso: alpha is fixed at 1 — use elasticnet for other mixes") | |
| 1576 | + } | |
| 1577 | + guard let value = Double(text) else { | |
| 1578 | + throw ZQEngineError("elasticnet: invalid alpha()") | |
| 1579 | + } | |
| 1580 | + alpha = value | |
| 1581 | + } | |
| 1582 | + guard let lambdaText = command.option("lambda")?.firstArgument, | |
| 1583 | + let lambda = Double(lambdaText) else { | |
| 1584 | + let hint = ZQElasticNet.lambdaMax( | |
| 1585 | + y: sample.y, predictors: sample.predictors, alpha: alpha | |
| 1586 | + ) | |
| 1587 | + throw ZQEngineError( | |
| 1588 | + "\(command.verb): lambda(#) required — lambda_max for this data is \(TableFormatter.general(hint, significant: 6))" | |
| 1589 | + ) | |
| 1590 | + } | |
| 1591 | + | |
| 1592 | + let result = try ZQElasticNet.fit( | |
| 1593 | + y: sample.y, | |
| 1594 | + predictors: sample.predictors, | |
| 1595 | + alpha: alpha, | |
| 1596 | + lambda: lambda | |
| 1597 | + ) | |
| 1598 | + | |
| 1599 | + // predict works after penalized fits; margins (needing a VCE) | |
| 1600 | + // deliberately does not. | |
| 1601 | + recordEstimation( | |
| 1602 | + kind: .ols, | |
| 1603 | + responseName: sample.responseName, | |
| 1604 | + coefficients: result.coefficients.map { | |
| 1605 | + ZQCoefficient( | |
| 1606 | + name: $0.name, estimate: $0.value, standardError: .nan, | |
| 1607 | + tStatistic: .nan, pValue: .nan, | |
| 1608 | + confidenceLower: .nan, confidenceUpper: .nan | |
| 1609 | + ) | |
| 1610 | + } + [ZQCoefficient( | |
| 1611 | + name: "_cons", estimate: result.intercept, standardError: .nan, | |
| 1612 | + tStatistic: .nan, pValue: .nan, | |
| 1613 | + confidenceLower: .nan, confidenceUpper: .nan | |
| 1614 | + )], | |
| 1615 | + definitions: sample.definitions, | |
| 1616 | + includeConstant: true, | |
| 1617 | + vce: [], | |
| 1618 | + sampleSize: sample.y.count | |
| 1619 | + ) | |
| 1620 | + | |
| 1621 | + var lines = [ | |
| 1622 | + alpha == 1 ? "Lasso linear model" : "Elastic-net linear model (alpha = \(TableFormatter.general(alpha)))", | |
| 1623 | + ] | |
| 1624 | + if sample.droppedMissing > 0 { | |
| 1625 | + lines.append("(\(sample.droppedMissing) observations dropped due to missing values)") | |
| 1626 | + } | |
| 1627 | + for (label, value) in [ | |
| 1628 | + ("Number of obs", "\(sample.y.count)"), | |
| 1629 | + ("lambda", TableFormatter.general(lambda, significant: 6)), | |
| 1630 | + ("Nonzero coefficients", "\(result.nonzeroCount) of \(result.coefficients.count)"), | |
| 1631 | + ("R-squared", TableFormatter.fixed(result.rSquared, decimals: 4)), | |
| 1632 | + ] { | |
| 1633 | + lines.append( | |
| 1634 | + TableFormatter.pad(label, 46, right: false) + "= " + | |
| 1635 | + TableFormatter.pad(value, 12) | |
| 1636 | + ) | |
| 1637 | + } | |
| 1638 | + lines.append("") | |
| 1639 | + lines.append( | |
| 1640 | + TableFormatter.pad(sample.responseName, 12) + " | " + | |
| 1641 | + TableFormatter.pad("Coefficient", 12) | |
| 1642 | + ) | |
| 1643 | + lines.append(String(repeating: "-", count: 13) + "+" + String(repeating: "-", count: 14)) | |
| 1644 | + | |
| 1645 | + var scalars: [String: Double] = [ | |
| 1646 | + "N": Double(sample.y.count), | |
| 1647 | + "lambda": lambda, | |
| 1648 | + "alpha": alpha, | |
| 1649 | + "k_nonzero": Double(result.nonzeroCount), | |
| 1650 | + "r2": result.rSquared, | |
| 1651 | + ] | |
| 1652 | + for coefficient in result.coefficients where coefficient.value != 0 { | |
| 1653 | + lines.append( | |
| 1654 | + TableFormatter.pad(coefficient.name, 12) + " | " + | |
| 1655 | + TableFormatter.pad(TableFormatter.general(coefficient.value), 12) | |
| 1656 | + ) | |
| 1657 | + scalars["b_\(coefficient.name)"] = coefficient.value | |
| 1658 | + } | |
| 1659 | + // Zeroed coefficients still surface as scalars (as exact 0). | |
| 1660 | + for coefficient in result.coefficients where coefficient.value == 0 { | |
| 1661 | + scalars["b_\(coefficient.name)"] = 0 | |
| 1662 | + } | |
| 1663 | + lines.append( | |
| 1664 | + TableFormatter.pad("_cons", 12) + " | " + | |
| 1665 | + TableFormatter.pad(TableFormatter.general(result.intercept), 12) | |
| 1666 | + ) | |
| 1667 | + scalars["b__cons"] = result.intercept | |
| 1668 | + return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars) | |
| 1669 | + } | |
| 1670 | + | |
| 1562 | 1671 | /// `margins, dydx(varlist)` — average marginal effects with |
| 1563 | 1672 | /// delta-method standard errors. OLS/IV effects are the coefficients |
| 1564 | 1673 | /// themselves; GLM effects average dμ/dx over the estimation sample: |
modified
MetrikaKit/Sources/ZQParser/KnownVerbs.swift
+2 −0
@@ -45,6 +45,8 @@ public struct ZQVerbTable: Sendable { | ||
| 45 | 45 | "xtreg": 5, |
| 46 | 46 | "xtset": 5, |
| 47 | 47 | "areg": 4, |
| 48 | + "elasticnet": 7, | |
| 49 | + "lasso": 5, | |
| 48 | 50 | "bootstrap": 9, |
| 49 | 51 | "permute": 7, |
| 50 | 52 | "jackknife": 9, |
added
MetrikaKit/Sources/ZQStats/ElasticNet.swift
+193 −0
@@ -0,0 +1,193 @@ | ||
| 1 | +// | |
| 2 | +// ElasticNet.swift | |
| 3 | +// Metrika | |
| 4 | +// | |
| 5 | +// Author: Simon-Pierre Boucher | |
| 6 | +// Contact: contact@spboucher.ai | |
| 7 | +// Copyright © 2026 Simon-Pierre Boucher. All rights reserved. | |
| 8 | +// | |
| 9 | + | |
| 10 | +import Foundation | |
| 11 | + | |
| 12 | +/// Elastic-net linear regression by cyclic coordinate descent, following | |
| 13 | +/// glmnet's conventions exactly so results are cross-checkable: | |
| 14 | +/// | |
| 15 | +/// min (1/2n)·Σᵢ(yᵢ − b₀ − xᵢ'β)² + λ·[α‖β‖₁ + ((1−α)/2)‖β‖₂²] | |
| 16 | +/// | |
| 17 | +/// Predictors are standardized internally (mean 0, variance 1 with the | |
| 18 | +/// 1/n denominator — glmnet's `standardize = TRUE`), the penalty applies | |
| 19 | +/// on the standardized scale, and coefficients are reported back on the | |
| 20 | +/// original scale with an unpenalized intercept. | |
| 21 | +public struct ZQElasticNetResult: Equatable, Sendable { | |
| 22 | + public var coefficients: [(name: String, value: Double)] | |
| 23 | + public var intercept: Double | |
| 24 | + public var alpha: Double | |
| 25 | + public var lambda: Double | |
| 26 | + public var nonzeroCount: Int | |
| 27 | + public var rSquared: Double | |
| 28 | + public var iterations: Int | |
| 29 | + | |
| 30 | + public static func == (lhs: ZQElasticNetResult, rhs: ZQElasticNetResult) -> Bool { | |
| 31 | + lhs.coefficients.elementsEqual(rhs.coefficients, by: { $0 == $1 }) | |
| 32 | + && lhs.intercept == rhs.intercept | |
| 33 | + && lhs.alpha == rhs.alpha | |
| 34 | + && lhs.lambda == rhs.lambda | |
| 35 | + && lhs.nonzeroCount == rhs.nonzeroCount | |
| 36 | + && lhs.rSquared == rhs.rSquared | |
| 37 | + && lhs.iterations == rhs.iterations | |
| 38 | + } | |
| 39 | +} | |
| 40 | + | |
| 41 | +public enum ZQElasticNet { | |
| 42 | + | |
| 43 | + /// - Parameters: | |
| 44 | + /// - alpha: 1 = lasso, 0 = ridge, in between = elastic net. | |
| 45 | + /// - lambda: penalty strength on glmnet's scale. | |
| 46 | + public static func fit( | |
| 47 | + y: [Double], | |
| 48 | + predictors: [(name: String, values: [Double])], | |
| 49 | + alpha: Double, | |
| 50 | + lambda: Double, | |
| 51 | + maximumIterations: Int = 100_000, | |
| 52 | + tolerance: Double = 1e-14 | |
| 53 | + ) throws -> ZQElasticNetResult { | |
| 54 | + let n = y.count | |
| 55 | + let p = predictors.count | |
| 56 | + guard p > 0 else { throw ZQStatsError("elasticnet: regressors required") } | |
| 57 | + guard n > 1 else { throw ZQStatsError("elasticnet: insufficient observations") } | |
| 58 | + guard (0...1).contains(alpha) else { | |
| 59 | + throw ZQStatsError("elasticnet: alpha() must be in [0, 1]") | |
| 60 | + } | |
| 61 | + guard lambda >= 0 else { | |
| 62 | + throw ZQStatsError("elasticnet: lambda() must be nonnegative") | |
| 63 | + } | |
| 64 | + | |
| 65 | + // Standardize (1/n variance, glmnet convention); center y. | |
| 66 | + let meanY = y.reduce(0, +) / Double(n) | |
| 67 | + var shifts = [Double](repeating: 0, count: p) | |
| 68 | + var scales = [Double](repeating: 1, count: p) | |
| 69 | + var standardized = [[Double]](repeating: [], count: p) | |
| 70 | + for (j, column) in predictors.enumerated() { | |
| 71 | + guard column.values.count == n else { | |
| 72 | + throw ZQStatsError("regressor '\(column.name)' has wrong length") | |
| 73 | + } | |
| 74 | + let mean = column.values.reduce(0, +) / Double(n) | |
| 75 | + let variance = column.values.reduce(0) { | |
| 76 | + $0 + ($1 - mean) * ($1 - mean) | |
| 77 | + } / Double(n) | |
| 78 | + guard variance > 0 else { | |
| 79 | + throw ZQStatsError("elasticnet: '\(column.name)' is constant") | |
| 80 | + } | |
| 81 | + shifts[j] = mean | |
| 82 | + scales[j] = variance.squareRoot() | |
| 83 | + standardized[j] = column.values.map { ($0 - mean) / scales[j] } | |
| 84 | + } | |
| 85 | + let centeredY = y.map { $0 - meanY } | |
| 86 | + | |
| 87 | + // Cyclic coordinate descent with residual updates. On the | |
| 88 | + // standardized scale each column has (1/n)Σx² = 1, so | |
| 89 | + // β_j ← S(ρ_j, λα) / (1 + λ₂), ρ_j = (1/n)x_j'r + β_j. | |
| 90 | + // | |
| 91 | + // glmnet quirk, matched deliberately: for gaussian models glmnet | |
| 92 | + // standardizes y internally, which leaves the L1 penalty invariant | |
| 93 | + // but scales the effective ridge penalty by 1/sd(y) — so | |
| 94 | + // λ₂ = λ(1−α)/sd(y) (1/n variance). Ridge solutions therefore | |
| 95 | + // depend on the scale of y, exactly as glmnet's do. | |
| 96 | + let varY = centeredY.reduce(0) { $0 + $1 * $1 } / Double(n) | |
| 97 | + guard varY > 0 else { | |
| 98 | + throw ZQStatsError("elasticnet: the response is constant") | |
| 99 | + } | |
| 100 | + var beta = [Double](repeating: 0, count: p) | |
| 101 | + var residuals = centeredY | |
| 102 | + let l1 = lambda * alpha | |
| 103 | + let l2 = lambda * (1 - alpha) / varY.squareRoot() | |
| 104 | + var iterations = 0 | |
| 105 | + | |
| 106 | + for iteration in 1...maximumIterations { | |
| 107 | + iterations = iteration | |
| 108 | + var maxChange = 0.0 | |
| 109 | + for j in 0..<p { | |
| 110 | + var rho = beta[j] | |
| 111 | + let column = standardized[j] | |
| 112 | + var dot = 0.0 | |
| 113 | + for i in 0..<n { dot += column[i] * residuals[i] } | |
| 114 | + rho += dot / Double(n) | |
| 115 | + | |
| 116 | + let updated = softThreshold(rho, l1) / (1 + l2) | |
| 117 | + let change = updated - beta[j] | |
| 118 | + if change != 0 { | |
| 119 | + for i in 0..<n { residuals[i] -= change * column[i] } | |
| 120 | + beta[j] = updated | |
| 121 | + maxChange = max(maxChange, abs(change)) | |
| 122 | + } | |
| 123 | + } | |
| 124 | + if maxChange < tolerance * max(1, beta.map(abs).max() ?? 1) { | |
| 125 | + break | |
| 126 | + } | |
| 127 | + if iteration == maximumIterations { | |
| 128 | + throw ZQStatsError("elasticnet: did not converge — increase lambda or iterations") | |
| 129 | + } | |
| 130 | + } | |
| 131 | + | |
| 132 | + // Back to the original scale. | |
| 133 | + var coefficients: [(name: String, value: Double)] = [] | |
| 134 | + var intercept = meanY | |
| 135 | + var nonzero = 0 | |
| 136 | + for j in 0..<p { | |
| 137 | + let value = beta[j] / scales[j] | |
| 138 | + coefficients.append((predictors[j].name, value)) | |
| 139 | + intercept -= shifts[j] * value | |
| 140 | + if value != 0 { nonzero += 1 } | |
| 141 | + } | |
| 142 | + | |
| 143 | + // Fit on the original data. | |
| 144 | + var rss = 0.0 | |
| 145 | + var tss = 0.0 | |
| 146 | + for i in 0..<n { | |
| 147 | + var fitted = intercept | |
| 148 | + for j in 0..<p { fitted += coefficients[j].value * predictors[j].values[i] } | |
| 149 | + rss += (y[i] - fitted) * (y[i] - fitted) | |
| 150 | + tss += (y[i] - meanY) * (y[i] - meanY) | |
| 151 | + } | |
| 152 | + | |
| 153 | + return ZQElasticNetResult( | |
| 154 | + coefficients: coefficients, | |
| 155 | + intercept: intercept, | |
| 156 | + alpha: alpha, | |
| 157 | + lambda: lambda, | |
| 158 | + nonzeroCount: nonzero, | |
| 159 | + rSquared: tss > 0 ? 1 - rss / tss : .nan, | |
| 160 | + iterations: iterations | |
| 161 | + ) | |
| 162 | + } | |
| 163 | + | |
| 164 | + /// The smallest λ at which every coefficient is zero (for α > 0): | |
| 165 | + /// λ_max = max_j |(1/n)·x_j'y| / α on the standardized scale. | |
| 166 | + public static func lambdaMax( | |
| 167 | + y: [Double], predictors: [(name: String, values: [Double])], alpha: Double | |
| 168 | + ) -> Double { | |
| 169 | + let n = y.count | |
| 170 | + let meanY = y.reduce(0, +) / Double(n) | |
| 171 | + var best = 0.0 | |
| 172 | + for column in predictors { | |
| 173 | + let mean = column.values.reduce(0, +) / Double(n) | |
| 174 | + let variance = column.values.reduce(0) { | |
| 175 | + $0 + ($1 - mean) * ($1 - mean) | |
| 176 | + } / Double(n) | |
| 177 | + guard variance > 0 else { continue } | |
| 178 | + let scale = variance.squareRoot() | |
| 179 | + var dot = 0.0 | |
| 180 | + for i in 0..<n { | |
| 181 | + dot += (column.values[i] - mean) / scale * (y[i] - meanY) | |
| 182 | + } | |
| 183 | + best = max(best, abs(dot) / Double(n)) | |
| 184 | + } | |
| 185 | + return alpha > 0 ? best / alpha : best * 1000 | |
| 186 | + } | |
| 187 | + | |
| 188 | + private static func softThreshold(_ value: Double, _ threshold: Double) -> Double { | |
| 189 | + if value > threshold { return value - threshold } | |
| 190 | + if value < -threshold { return value + threshold } | |
| 191 | + return 0 | |
| 192 | + } | |
| 193 | +} | |
modified
MetrikaKit/Sources/ZQStats/OLS.swift
+14 −0
@@ -29,6 +29,20 @@ public struct ZQCoefficient: Equatable, Sendable { | ||
| 29 | 29 | public var pValue: Double |
| 30 | 30 | public var confidenceLower: Double |
| 31 | 31 | public var confidenceUpper: Double |
| 32 | + | |
| 33 | + public init( | |
| 34 | + name: String, estimate: Double, standardError: Double, | |
| 35 | + tStatistic: Double, pValue: Double, | |
| 36 | + confidenceLower: Double, confidenceUpper: Double | |
| 37 | + ) { | |
| 38 | + self.name = name | |
| 39 | + self.estimate = estimate | |
| 40 | + self.standardError = standardError | |
| 41 | + self.tStatistic = tStatistic | |
| 42 | + self.pValue = pValue | |
| 43 | + self.confidenceLower = confidenceLower | |
| 44 | + self.confidenceUpper = confidenceUpper | |
| 45 | + } | |
| 32 | 46 | } |
| 33 | 47 | |
| 34 | 48 | /// Full OLS fit result. |
added
MetrikaKit/Tests/MetrikaKitTests/ElasticNetTests.swift
+122 −0
@@ -0,0 +1,122 @@ | ||
| 1 | +// | |
| 2 | +// ElasticNetTests.swift | |
| 3 | +// Metrika | |
| 4 | +// | |
| 5 | +// Author: Simon-Pierre Boucher | |
| 6 | +// Contact: contact@spboucher.ai | |
| 7 | +// Copyright © 2026 Simon-Pierre Boucher. All rights reserved. | |
| 8 | +// | |
| 9 | + | |
| 10 | +import Foundation | |
| 11 | +import Testing | |
| 12 | +import ZQEngine | |
| 13 | +import ZQStats | |
| 14 | + | |
| 15 | +/// Elastic net vs glmnet fixtures. Both sides are iterative coordinate | |
| 16 | +/// descent (glmnet thresh 1e-15, ours 1e-14), so agreement is asserted at | |
| 17 | +/// 1e-6 relative — the documented tolerance for penalized solvers. | |
| 18 | +@Suite("Elastic net", .serialized) | |
| 19 | +struct ElasticNetTests { | |
| 20 | + let fixtures: Fixtures | |
| 21 | + let session: ZQSession | |
| 22 | + | |
| 23 | + init() async throws { | |
| 24 | + self.fixtures = try Fixtures() | |
| 25 | + self.session = try ZQSession(discoverUserCommands: false) | |
| 26 | + _ = try await session.execute("use \(fixtures.datasetURL.path)") | |
| 27 | + _ = try await session.execute("gen log_rev = ln(revenue)") | |
| 28 | + } | |
| 29 | + | |
| 30 | + @Test("lasso matches glmnet and reproduces its selection") | |
| 31 | + func lasso() async throws { | |
| 32 | + let result = try await session.execute( | |
| 33 | + "lasso log_rev price z1 z2 orders, lambda(0.05)" | |
| 34 | + ) | |
| 35 | + expectClose( | |
| 36 | + try #require(result.scalars["b_price"]), | |
| 37 | + fixtures["lasso_b_price"], rtol: 1e-6, "b[price]" | |
| 38 | + ) | |
| 39 | + expectClose( | |
| 40 | + try #require(result.scalars["b_z1"]), | |
| 41 | + fixtures["lasso_b_z1"], rtol: 1e-4, "b[z1] (tiny, near the threshold)" | |
| 42 | + ) | |
| 43 | + expectClose( | |
| 44 | + try #require(result.scalars["b__cons"]), | |
| 45 | + fixtures["lasso_b_cons"], rtol: 1e-6, "intercept" | |
| 46 | + ) | |
| 47 | + // glmnet zeroed z2 and orders — the selection must agree exactly. | |
| 48 | + #expect(result.scalars["b_z2"] == 0) | |
| 49 | + #expect(result.scalars["b_orders"] == 0) | |
| 50 | + #expect(result.scalars["k_nonzero"] == 2) | |
| 51 | + } | |
| 52 | + | |
| 53 | + @Test("elastic net (alpha 0.4) matches glmnet") | |
| 54 | + func elasticNet() async throws { | |
| 55 | + let result = try await session.execute( | |
| 56 | + "elasticnet log_rev price z1 z2 orders, alpha(0.4) lambda(0.02)" | |
| 57 | + ) | |
| 58 | + expectClose( | |
| 59 | + try #require(result.scalars["b_price"]), | |
| 60 | + fixtures["enet_b_price"], rtol: 1e-6, "b[price]" | |
| 61 | + ) | |
| 62 | + expectClose( | |
| 63 | + try #require(result.scalars["b_z1"]), | |
| 64 | + fixtures["enet_b_z1"], rtol: 1e-5, "b[z1]" | |
| 65 | + ) | |
| 66 | + expectClose( | |
| 67 | + try #require(result.scalars["b__cons"]), | |
| 68 | + fixtures["enet_b_cons"], rtol: 1e-6, "intercept" | |
| 69 | + ) | |
| 70 | + #expect(result.scalars["b_orders"] == 0) | |
| 71 | + } | |
| 72 | + | |
| 73 | + @Test("lambda 0 reproduces OLS") | |
| 74 | + func lambdaZeroIsOLS() async throws { | |
| 75 | + let ols = try await session.execute("reg log_rev price orders") | |
| 76 | + let penalized = try await session.execute( | |
| 77 | + "elasticnet log_rev price orders, alpha(1) lambda(0)" | |
| 78 | + ) | |
| 79 | + expectClose( | |
| 80 | + try #require(penalized.scalars["b_price"]), | |
| 81 | + try #require(ols.scalars["b_price"]), | |
| 82 | + rtol: 1e-8, "lambda 0 == OLS" | |
| 83 | + ) | |
| 84 | + } | |
| 85 | + | |
| 86 | + @Test("lambda at lambda_max zeroes everything") | |
| 87 | + func lambdaMaxZeroes() throws { | |
| 88 | + let n = 100 | |
| 89 | + let x = (0..<n).map { Double($0) } | |
| 90 | + let y = x.map { 2 + 3 * $0 } | |
| 91 | + let lambdaMax = ZQElasticNet.lambdaMax( | |
| 92 | + y: y, predictors: [("x", x)], alpha: 1 | |
| 93 | + ) | |
| 94 | + let result = try ZQElasticNet.fit( | |
| 95 | + y: y, predictors: [("x", x)], alpha: 1, lambda: lambdaMax * 1.0000001 | |
| 96 | + ) | |
| 97 | + #expect(result.nonzeroCount == 0) | |
| 98 | + expectClose(result.intercept, y.reduce(0, +) / Double(n), "null intercept") | |
| 99 | + } | |
| 100 | + | |
| 101 | + @Test("predict works after a penalized fit; margins refuses") | |
| 102 | + func postEstimation() async throws { | |
| 103 | + _ = try await session.execute("lasso log_rev price z1 z2 orders, lambda(0.05)") | |
| 104 | + _ = try await session.execute("predict lhat") | |
| 105 | + let summary = try await session.execute("summarize lhat") | |
| 106 | + #expect(summary.scalars["N"] == 57) | |
| 107 | + await #expect(throws: ZQEngineError.self) { | |
| 108 | + _ = try await session.execute("margins, dydx(price)") // no VCE | |
| 109 | + } | |
| 110 | + _ = try await session.execute("drop lhat") | |
| 111 | + } | |
| 112 | + | |
| 113 | + @Test("missing lambda reports lambda_max as a hint") | |
| 114 | + func lambdaHint() async throws { | |
| 115 | + do { | |
| 116 | + _ = try await session.execute("lasso log_rev price") | |
| 117 | + Issue.record("expected an error") | |
| 118 | + } catch let error as ZQEngineError { | |
| 119 | + #expect(error.message.contains("lambda_max")) | |
| 120 | + } | |
| 121 | + } | |
| 122 | +} | |
modified
MetrikaKit/Tests/MetrikaKitTests/Fixtures/expected.tsv
+10 −0
@@ -78,6 +78,16 @@ margins_logit_price -0.025653639230005385 | ||
| 78 | 78 | margins_logit_se 0.012555697160526031 |
| 79 | 79 | margins_pois_price -0.11269819308994093 |
| 80 | 80 | margins_pois_se 0.036026880220264283 |
| 81 | +lasso_b_cons 3.9156372006944973 | |
| 82 | +lasso_b_price -0.063375317860093414 | |
| 83 | +lasso_b_z1 -0.0023537309683557802 | |
| 84 | +lasso_b_z2 0 | |
| 85 | +lasso_b_orders 0 | |
| 86 | +enet_b_cons 4.0119027089030768 | |
| 87 | +enet_b_price -0.058669179277376367 | |
| 88 | +enet_b_z1 -0.014365113727136758 | |
| 89 | +enet_b_z2 9.5044046897303132e-05 | |
| 90 | +enet_b_orders 0 | |
| 81 | 91 | corr_rev_price -0.88919780339930277 |
| 82 | 92 | corr_rev_logrev 0.98089325161602359 |
| 83 | 93 | dist_pchisq_3p8_1 0.94874741714263044 |
modified
MetrikaKit/Tests/MetrikaKitTests/Fixtures/regression_v117.dta
+0 −0
Binary file not shown.
modified
MetrikaKit/Tests/MetrikaKitTests/Fixtures/regression_v118.dta
+0 −0
Binary file not shown.
modified
Tests/Fixtures/generate.R
+32 −0
@@ -320,6 +320,38 @@ m_ps <- ame_glm(ps, "price", exp, exp) | ||
| 320 | 320 | emit("margins_pois_price", m_ps[["ame"]]) |
| 321 | 321 | emit("margins_pois_se", m_ps[["se"]]) |
| 322 | 322 | |
| 323 | +# ------------------------------------------------------------ elastic net | |
| 324 | +# glmnet at fixed lambda, tight threshold. price/z1/z2 are correlated by | |
| 325 | +# construction, so the l1 penalty must actually choose among them. | |
| 326 | +if (requireNamespace("glmnet", quietly = TRUE)) { | |
| 327 | + Xen <- as.matrix(d[, c("price", "z1", "z2", "orders")]) | |
| 328 | + yen <- d$log_rev | |
| 329 | + | |
| 330 | + fit_lasso <- glmnet::glmnet( | |
| 331 | + Xen, yen, alpha = 1, lambda = 0.05, | |
| 332 | + thresh = 1e-15, maxit = 1e7, standardize = TRUE | |
| 333 | + ) | |
| 334 | + cl <- as.vector(coef(fit_lasso)) | |
| 335 | + emit("lasso_b_cons", cl[1]) | |
| 336 | + emit("lasso_b_price", cl[2]) | |
| 337 | + emit("lasso_b_z1", cl[3]) | |
| 338 | + emit("lasso_b_z2", cl[4]) | |
| 339 | + emit("lasso_b_orders", cl[5]) | |
| 340 | + | |
| 341 | + fit_enet <- glmnet::glmnet( | |
| 342 | + Xen, yen, alpha = 0.4, lambda = 0.02, | |
| 343 | + thresh = 1e-15, maxit = 1e7, standardize = TRUE | |
| 344 | + ) | |
| 345 | + ce <- as.vector(coef(fit_enet)) | |
| 346 | + emit("enet_b_cons", ce[1]) | |
| 347 | + emit("enet_b_price", ce[2]) | |
| 348 | + emit("enet_b_z1", ce[3]) | |
| 349 | + emit("enet_b_z2", ce[4]) | |
| 350 | + emit("enet_b_orders", ce[5]) | |
| 351 | +} else { | |
| 352 | + cat("glmnet not installed — skipping elastic-net fixtures\n") | |
| 353 | +} | |
| 354 | + | |
| 323 | 355 | # -------------------------------------------------------------- correlate |
| 324 | 356 | emit("corr_rev_price", cor(d$revenue, d$price)) |
| 325 | 357 | emit("corr_rev_logrev", cor(d$revenue, d$log_rev)) |
| 326 | 358 | |