SPB Git
spboucher.ai

Repository Public

spbmetrika

Commit

feat(engine): margins dydx() with delta-method standard errors

- ZQOLSResult/ZQGLMResult expose the full covariance matrix (column-major
  k x k aligned with coefficients); IV passes it through
- EstimationState carries vce, inference df, and the estimation-sample
  design (GLMs only — AMEs need it)
- margins, dydx(varlist): OLS/IV effects are the coefficients with their
  SEs; GLM average marginal effects with analytic delta gradients
  (logit p(1-p)(1-2p), probit -xb*phi, poisson exp) over the estimation
  sample; t or z inference per model kind
- factor/interaction dydx rejected with a clear message (discrete-change
  margins later); continuous terms of factor models work
- R fixtures mirror the exact formulas at the converged coefficients;
  logit/poisson AME and delta SE match at 1e-10
- 98 tests green (swift test and xcodebuild with GPU suites)

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
simon-pierre boucher committed · Aug 5, 2026 parent 1f3fd3b
10 changed files +349 −10 Browse files
10 changed files
Unified Split
modified MetrikaKit/Sources/ZQEngine/Session.swift +208 −5
@@ -195,6 +195,7 @@ public actor ZQSession {
195 case "xtreg": return try handleXTReg(command) 195 case "xtreg": return try handleXTReg(command)
196 case "ivregress": return try handleIVRegress(command) 196 case "ivregress": return try handleIVRegress(command)
197 case "predict": return try handlePredict(command) 197 case "predict": return try handlePredict(command)
198+ case "margins": return try handleMargins(command)
198 case "logit": return try handleGLM(command, family: .logit) 199 case "logit": return try handleGLM(command, family: .logit)
199 case "probit": return try handleGLM(command, family: .probit) 200 case "probit": return try handleGLM(command, family: .probit)
200 case "poisson": return try handleGLM(command, family: .poisson) 201 case "poisson": return try handleGLM(command, family: .poisson)
@@ -554,6 +555,14 @@ public actor ZQSession {
554 var responseName: String 555 var responseName: String
555 /// Coefficients aligned with their recomputation recipes. 556 /// Coefficients aligned with their recomputation recipes.
556 var coefficients: [(name: String, value: Double, definition: PredictorDefinition)] 557 var coefficients: [(name: String, value: Double, definition: PredictorDefinition)]
558+ /// Covariance matrix, column-major k×k aligned with `coefficients`.
559+ var vce: [Double]
560+ /// Degrees of freedom for t inference; nil means normal (ML).
561+ var inferenceDF: Double?
562+ /// Estimation-sample design (column-major n×k incl. constant) —
563+ /// kept for GLMs, whose average marginal effects need it.
564+ var sampleDesign: [Double]?
565+ var sampleSize: Int
557 } 566 }
558 567
559 /// Assembled estimation sample after if/in restriction and listwise 568 /// Assembled estimation sample after if/in restriction and listwise
@@ -740,7 +749,11 @@ public actor ZQSession {
740 responseName: String, 749 responseName: String,
741 coefficients: [ZQCoefficient], 750 coefficients: [ZQCoefficient],
742 definitions: [PredictorDefinition], 751 definitions: [PredictorDefinition],
743− includeConstant: Bool 752+ includeConstant: Bool,
753+ vce: [Double] = [],
754+ inferenceDF: Double? = nil,
755+ sampleDesign: [Double]? = nil,
756+ sampleSize: Int = 0
744 ) { 757 ) {
745 var recipes = definitions 758 var recipes = definitions
746 if includeConstant { recipes.append(.constant) } 759 if includeConstant { recipes.append(.constant) }
@@ -753,10 +766,26 @@ public actor ZQSession {
753 responseName: responseName, 766 responseName: responseName,
754 coefficients: zip(coefficients, recipes).map { 767 coefficients: zip(coefficients, recipes).map {
755 (name: $0.name, value: $0.estimate, definition: $1) 768 (name: $0.name, value: $0.estimate, definition: $1)
756− } 769+ },
770+ vce: vce,
771+ inferenceDF: inferenceDF,
772+ sampleDesign: sampleDesign,
773+ sampleSize: sampleSize
757 ) 774 )
758 } 775 }
759 776
777+ /// Column-major estimation-sample design from sample predictors.
778+ private func designMatrix(
779+ from sample: RegressionSample, includeConstant: Bool
780+ ) -> [Double] {
781+ let n = sample.y.count
782+ var design = [Double]()
783+ design.reserveCapacity(n * (sample.predictors.count + 1))
784+ for column in sample.predictors { design.append(contentsOf: column.values) }
785+ if includeConstant { design.append(contentsOf: [Double](repeating: 1, count: n)) }
786+ return design
787+ }
788+
760 private func varianceEstimator( 789 private func varianceEstimator(
761 _ command: ZQCommand, clusterLabels: [Int]? 790 _ command: ZQCommand, clusterLabels: [Int]?
762 ) throws -> ZQVarianceEstimator { 791 ) throws -> ZQVarianceEstimator {
@@ -795,7 +824,10 @@ public actor ZQSession {
795 responseName: sample.responseName, 824 responseName: sample.responseName,
796 coefficients: result.coefficients, 825 coefficients: result.coefficients,
797 definitions: sample.definitions, 826 definitions: sample.definitions,
798− includeConstant: !command.hasOption("noconstant") 827+ includeConstant: !command.hasOption("noconstant"),
828+ vce: result.vce,
829+ inferenceDF: result.inferenceDF,
830+ sampleSize: result.observationCount
799 ) 831 )
800 832
801 var header: [String] = [] 833 var header: [String] = []
@@ -1023,7 +1055,10 @@ public actor ZQSession {
1023 responseName: response, 1055 responseName: response,
1024 coefficients: result.coefficients, 1056 coefficients: result.coefficients,
1025 definitions: (ivSpec.endogenous + exogenousNames).map { .column($0) }, 1057 definitions: (ivSpec.endogenous + exogenousNames).map { .column($0) },
1026− includeConstant: !command.hasOption("noconstant") 1058+ includeConstant: !command.hasOption("noconstant"),
1059+ vce: result.vce,
1060+ inferenceDF: result.inferenceDF,
1061+ sampleSize: result.observationCount
1027 ) 1062 )
1028 1063
1029 var header = ["Instrumental variables (2SLS) regression"] 1064 var header = ["Instrumental variables (2SLS) regression"]
@@ -1125,7 +1160,13 @@ public actor ZQSession {
1125 responseName: sample.responseName, 1160 responseName: sample.responseName,
1126 coefficients: result.coefficients, 1161 coefficients: result.coefficients,
1127 definitions: sample.definitions, 1162 definitions: sample.definitions,
1128− includeConstant: !command.hasOption("noconstant") 1163+ includeConstant: !command.hasOption("noconstant"),
1164+ vce: result.vce,
1165+ sampleDesign: designMatrix(
1166+ from: sample,
1167+ includeConstant: !command.hasOption("noconstant")
1168+ ),
1169+ sampleSize: result.observationCount
1129 ) 1170 )
1130 1171
1131 let title: String 1172 let title: String
@@ -1518,6 +1559,168 @@ public actor ZQSession {
1518 return ZQResult(text: note, scalars: ["N": Double(n - missingCount)]) 1559 return ZQResult(text: note, scalars: ["N": Double(n - missingCount)])
1519 } 1560 }
1520 1561
1562+ /// `margins, dydx(varlist)` — average marginal effects with
1563+ /// delta-method standard errors. OLS/IV effects are the coefficients
1564+ /// themselves; GLM effects average dμ/dx over the estimation sample:
1565+ /// AME_j = β_j · (1/N) Σᵢ g′(xbᵢ)
1566+ /// ∂AME_j/∂β_m = (1/N) Σᵢ g″(xbᵢ)·x_im·β_j + δ_jm·(1/N) Σᵢ g′(xbᵢ)
1567+ /// Only continuous (plain-column) regressors are supported so far.
1568+ private func handleMargins(_ command: ZQCommand) throws -> ZQResult {
1569+ guard let estimation = lastEstimation else {
1570+ throw ZQEngineError("margins: no estimation results — run a regression first")
1571+ }
1572+ guard let dydx = command.option("dydx"), !dydx.arguments.isEmpty else {
1573+ throw ZQEngineError("margins: syntax is 'margins, dydx(varlist)'")
1574+ }
1575+ guard !estimation.vce.isEmpty else {
1576+ throw ZQEngineError("margins: the last estimation stored no covariance matrix")
1577+ }
1578+
1579+ let k = estimation.coefficients.count
1580+ struct Effect {
1581+ var name: String
1582+ var dydx: Double
1583+ var standardError: Double
1584+ }
1585+ var effects: [Effect] = []
1586+
1587+ for name in dydx.arguments {
1588+ guard let j = estimation.coefficients.firstIndex(where: {
1589+ $0.definition == .column(name)
1590+ }) else {
1591+ if estimation.coefficients.contains(where: {
1592+ if case .indicator(let variable, _) = $0.definition { return variable == name }
1593+ if case .product(let variables) = $0.definition { return variables.contains(name) }
1594+ return false
1595+ }) {
1596+ throw ZQEngineError(
1597+ "margins: dydx over factor or interaction terms is not supported yet"
1598+ )
1599+ }
1600+ throw ZQEngineError("margins: '\(name)' is not a regressor in the last estimation")
1601+ }
1602+ let beta = estimation.coefficients[j].value
1603+
1604+ switch estimation.kind {
1605+ case .ols, .iv:
1606+ // Linear model: the marginal effect is the coefficient.
1607+ effects.append(Effect(
1608+ name: name,
1609+ dydx: beta,
1610+ standardError: estimation.vce[j * k + j].squareRoot()
1611+ ))
1612+
1613+ case .glm(let family):
1614+ guard let design = estimation.sampleDesign,
1615+ estimation.sampleSize > 0 else {
1616+ throw ZQEngineError("margins: estimation sample unavailable")
1617+ }
1618+ let n = estimation.sampleSize
1619+ // xb over the estimation sample.
1620+ var xb = [Double](repeating: 0, count: n)
1621+ for m in 0..<k {
1622+ let value = estimation.coefficients[m].value
1623+ for i in 0..<n { xb[i] += value * design[m * n + i] }
1624+ }
1625+ // g′ and g″ per observation.
1626+ var meanFirst = 0.0
1627+ var second = [Double](repeating: 0, count: n)
1628+ for i in 0..<n {
1629+ switch family {
1630+ case .logit:
1631+ let p = 1 / (1 + Foundation.exp(-xb[i]))
1632+ meanFirst += p * (1 - p)
1633+ second[i] = p * (1 - p) * (1 - 2 * p)
1634+ case .probit:
1635+ let density = Foundation.exp(-0.5 * xb[i] * xb[i])
1636+ / (2 * Double.pi).squareRoot()
1637+ meanFirst += density
1638+ second[i] = -xb[i] * density
1639+ case .poisson:
1640+ let mu = Foundation.exp(xb[i])
1641+ meanFirst += mu
1642+ second[i] = mu
1643+ }
1644+ }
1645+ meanFirst /= Double(n)
1646+
1647+ // Delta-method gradient.
1648+ var gradient = [Double](repeating: 0, count: k)
1649+ for m in 0..<k {
1650+ var meanSecondX = 0.0
1651+ for i in 0..<n { meanSecondX += second[i] * design[m * n + i] }
1652+ gradient[m] = beta * meanSecondX / Double(n)
1653+ }
1654+ gradient[j] += meanFirst
1655+
1656+ var variance = 0.0
1657+ for a in 0..<k {
1658+ for b in 0..<k {
1659+ variance += gradient[a] * estimation.vce[b * k + a] * gradient[b]
1660+ }
1661+ }
1662+ effects.append(Effect(
1663+ name: name,
1664+ dydx: beta * meanFirst,
1665+ standardError: variance.squareRoot()
1666+ ))
1667+ }
1668+ }
1669+
1670+ // Render: t inference when the model carries a df, z otherwise.
1671+ let level = try confidenceLevel(command)
1672+ let usesT = estimation.inferenceDF != nil
1673+ let critical: Double
1674+ if let df = estimation.inferenceDF {
1675+ critical = ZQDistributions.studentTQuantile(0.5 + level / 2, df: df)
1676+ } else {
1677+ critical = ZQDistributions.normalQuantile(0.5 + level / 2)
1678+ }
1679+
1680+ var lines = [
1681+ "Average marginal effects Number of obs = " +
1682+ TableFormatter.pad("\(estimation.sampleSize)", 10),
1683+ "",
1684+ ]
1685+ let widths = [12, 12, 11, 8, 8, 22]
1686+ let statLabel = usesT ? "t" : "z"
1687+ lines.append(
1688+ TableFormatter.pad("", widths[0]) + " | " +
1689+ TableFormatter.pad("dy/dx", widths[1]) + " " +
1690+ TableFormatter.pad("Std. err.", widths[2]) + " " +
1691+ TableFormatter.pad(statLabel, widths[3]) + " " +
1692+ TableFormatter.pad("P>|\(statLabel)|", widths[4]) + " " +
1693+ TableFormatter.pad("[\(Int(level * 100))% conf. interval]", widths[5])
1694+ )
1695+ lines.append(TableFormatter.rule(widths))
1696+
1697+ var scalars: [String: Double] = ["N": Double(estimation.sampleSize)]
1698+ for effect in effects {
1699+ let statistic = effect.dydx / effect.standardError
1700+ let p: Double
1701+ if let df = estimation.inferenceDF {
1702+ p = ZQDistributions.tTestPValue(statistic, df: df)
1703+ } else {
1704+ p = 2 * (1 - ZQDistributions.normalCDF(abs(statistic)))
1705+ }
1706+ lines.append(
1707+ TableFormatter.pad(effect.name, widths[0]) + " | " +
1708+ TableFormatter.pad(TableFormatter.general(effect.dydx), widths[1]) + " " +
1709+ TableFormatter.pad(TableFormatter.general(effect.standardError), widths[2]) + " " +
1710+ TableFormatter.pad(TableFormatter.fixed(statistic, decimals: 2), widths[3]) + " " +
1711+ TableFormatter.pad(TableFormatter.fixed(p, decimals: 3), widths[4]) + " " +
1712+ TableFormatter.pad(
1713+ TableFormatter.general(effect.dydx - critical * effect.standardError) + " " +
1714+ TableFormatter.general(effect.dydx + critical * effect.standardError),
1715+ widths[5]
1716+ )
1717+ )
1718+ scalars["dydx_\(effect.name)"] = effect.dydx
1719+ scalars["se_\(effect.name)"] = effect.standardError
1720+ }
1721+ return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
1722+ }
1723+
1521 // MARK: - Bootstrap prefix 1724 // MARK: - Bootstrap prefix
1522 1725
1523 private func handleBootstrap( 1726 private func handleBootstrap(
modified MetrikaKit/Sources/ZQParser/KnownVerbs.swift +1 −0
@@ -32,6 +32,7 @@ public struct ZQVerbTable: Sendable {
32 "correlate": 3, 32 "correlate": 3,
33 "describe": 1, 33 "describe": 1,
34 "list": 1, 34 "list": 1,
35+ "margins": 4,
35 "generate": 3, 36 "generate": 3,
36 "replace": 7, 37 "replace": 7,
37 "drop": 4, 38 "drop": 4,
modified MetrikaKit/Sources/ZQStats/GLM.swift +8 −2
@@ -34,6 +34,9 @@ public struct ZQGLMResult: Equatable, Sendable {
34 public var clusterCount: Int? 34 public var clusterCount: Int?
35 /// Linear predictor and fitted mean at the optimum. 35 /// Linear predictor and fitted mean at the optimum.
36 public var fittedMeans: [Double] 36 public var fittedMeans: [Double]
37+ /// Full covariance matrix, column-major k×k, aligned with
38+ /// `coefficients`.
39+ public var vce: [Double]
37 40
38 public init( 41 public init(
39 family: ZQGLMFamily, 42 family: ZQGLMFamily,
@@ -47,7 +50,8 @@ public struct ZQGLMResult: Equatable, Sendable {
47 chiSquaredPValue: Double?, 50 chiSquaredPValue: Double?,
48 iterations: Int, 51 iterations: Int,
49 clusterCount: Int?, 52 clusterCount: Int?,
50− fittedMeans: [Double] 53+ fittedMeans: [Double],
54+ vce: [Double] = []
51 ) { 55 ) {
52 self.family = family 56 self.family = family
53 self.coefficients = coefficients 57 self.coefficients = coefficients
@@ -61,6 +65,7 @@ public struct ZQGLMResult: Equatable, Sendable {
61 self.iterations = iterations 65 self.iterations = iterations
62 self.clusterCount = clusterCount 66 self.clusterCount = clusterCount
63 self.fittedMeans = fittedMeans 67 self.fittedMeans = fittedMeans
68+ self.vce = vce
64 } 69 }
65 } 70 }
66 71
@@ -324,7 +329,8 @@ public enum ZQGLM {
324 chiSquaredPValue: chiSquaredPValue, 329 chiSquaredPValue: chiSquaredPValue,
325 iterations: iterations, 330 iterations: iterations,
326 clusterCount: clusterCount, 331 clusterCount: clusterCount,
327− fittedMeans: mu 332+ fittedMeans: mu,
333+ vce: vce
328 ) 334 )
329 } 335 }
330 336
modified MetrikaKit/Sources/ZQStats/IV.swift +2 −1
@@ -228,7 +228,8 @@ public enum ZQIV {
228 fPValue: fPValue, 228 fPValue: fPValue,
229 fDF: fDF, 229 fDF: fDF,
230 clusterCount: clusterCount, 230 clusterCount: clusterCount,
231− residuals: residuals 231+ residuals: residuals,
232+ vce: vce
232 ) 233 )
233 } 234 }
234 } 235 }
modified MetrikaKit/Sources/ZQStats/OLS.swift +9 −2
@@ -50,6 +50,10 @@ public struct ZQOLSResult: Equatable, Sendable {
50 } 50 }
51 public var clusterCount: Int? 51 public var clusterCount: Int?
52 public var residuals: [Double] 52 public var residuals: [Double]
53+ /// Full covariance matrix of the coefficients, column-major k×k,
54+ /// aligned with `coefficients` — needed by delta-method
55+ /// post-estimation (margins).
56+ public var vce: [Double]
53 57
54 private struct FDF: Equatable, Sendable { 58 private struct FDF: Equatable, Sendable {
55 let first: Double 59 let first: Double
@@ -75,7 +79,8 @@ public struct ZQOLSResult: Equatable, Sendable {
75 fPValue: Double?, 79 fPValue: Double?,
76 fDF: (Double, Double)?, 80 fDF: (Double, Double)?,
77 clusterCount: Int?, 81 clusterCount: Int?,
78− residuals: [Double] 82+ residuals: [Double],
83+ vce: [Double] = []
79 ) { 84 ) {
80 self.coefficients = coefficients 85 self.coefficients = coefficients
81 self.observationCount = observationCount 86 self.observationCount = observationCount
@@ -89,6 +94,7 @@ public struct ZQOLSResult: Equatable, Sendable {
89 self.fDFStorage = fDF.map { FDF($0.0, $0.1) } 94 self.fDFStorage = fDF.map { FDF($0.0, $0.1) }
90 self.clusterCount = clusterCount 95 self.clusterCount = clusterCount
91 self.residuals = residuals 96 self.residuals = residuals
97+ self.vce = vce
92 } 98 }
93 } 99 }
94 100
@@ -292,7 +298,8 @@ public enum ZQOLS {
292 fPValue: fPValue, 298 fPValue: fPValue,
293 fDF: fDF, 299 fDF: fDF,
294 clusterCount: clusterCount, 300 clusterCount: clusterCount,
295− residuals: residuals 301+ residuals: residuals,
302+ vce: vce
296 ) 303 )
297 } 304 }
298 305
modified MetrikaKit/Tests/MetrikaKitTests/Fixtures/expected.tsv +4 −0
@@ -74,6 +74,10 @@ iv_p_price 5.7814855186358344e-21
74 iv_r2 0.81943731431549249 74 iv_r2 0.81943731431549249
75 iv_rmse 0.16468550892743888 75 iv_rmse 0.16468550892743888
76 iv_se_hc1_price 0.005144683256992212 76 iv_se_hc1_price 0.005144683256992212
77+margins_logit_price -0.025653639230005385
78+margins_logit_se 0.012555697160526031
79+margins_pois_price -0.11269819308994093
80+margins_pois_se 0.036026880220264283
77 corr_rev_price -0.88919780339930277 81 corr_rev_price -0.88919780339930277
78 corr_rev_logrev 0.98089325161602359 82 corr_rev_logrev 0.98089325161602359
79 dist_pchisq_3p8_1 0.94874741714263044 83 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.

added MetrikaKit/Tests/MetrikaKitTests/MarginsTests.swift +88 −0
@@ -0,0 +1,88 @@
1+//
2+// MarginsTests.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+
14+@Suite("margins", .serialized)
15+struct MarginsTests {
16+ let fixtures: Fixtures
17+ let session: ZQSession
18+
19+ init() async throws {
20+ self.fixtures = try Fixtures()
21+ self.session = try ZQSession(discoverUserCommands: false)
22+ _ = try await session.execute("use \(fixtures.datasetURL.path)")
23+ _ = try await session.execute("gen log_rev = ln(revenue)")
24+ }
25+
26+ @Test("OLS marginal effect is the coefficient with its SE")
27+ func olsMargins() async throws {
28+ let fit = try await session.execute("reg log_rev price")
29+ let margins = try await session.execute("margins, dydx(price)")
30+ #expect(margins.scalars["dydx_price"] == fit.scalars["b_price"])
31+ #expect(margins.scalars["se_price"] == fit.scalars["se_price"])
32+ }
33+
34+ @Test("logit AME and delta-method SE match R")
35+ func logitMargins() async throws {
36+ _ = try await session.execute("logit purchase price")
37+ let margins = try await session.execute("margins, dydx(price)")
38+ expectClose(
39+ try #require(margins.scalars["dydx_price"]),
40+ fixtures["margins_logit_price"], "logit AME"
41+ )
42+ expectClose(
43+ try #require(margins.scalars["se_price"]),
44+ fixtures["margins_logit_se"], "logit delta SE"
45+ )
46+ }
47+
48+ @Test("poisson AME and delta-method SE match R")
49+ func poissonMargins() async throws {
50+ _ = try await session.execute("poisson orders price")
51+ let margins = try await session.execute("margins, dydx(price)")
52+ expectClose(
53+ try #require(margins.scalars["dydx_price"]),
54+ fixtures["margins_pois_price"], "poisson AME"
55+ )
56+ expectClose(
57+ try #require(margins.scalars["se_price"]),
58+ fixtures["margins_pois_se"], "poisson delta SE"
59+ )
60+ }
61+
62+ @Test("dydx over factor terms is rejected for now")
63+ func factorRejected() async throws {
64+ _ = try await session.execute("reg log_rev price i.region")
65+ await #expect(throws: ZQEngineError.self) {
66+ _ = try await session.execute("margins, dydx(region)")
67+ }
68+ // Continuous term of the same model still works.
69+ let margins = try await session.execute("margins, dydx(price)")
70+ #expect(margins.scalars["dydx_price"] != nil)
71+ }
72+
73+ @Test("margins requires estimation results and a dydx() option")
74+ func validation() async throws {
75+ let fresh = try ZQSession(discoverUserCommands: false)
76+ _ = try await fresh.execute("use \(fixtures.datasetURL.path)")
77+ await #expect(throws: ZQEngineError.self) {
78+ _ = try await fresh.execute("margins, dydx(price)")
79+ }
80+ _ = try await session.execute("reg log_rev price")
81+ await #expect(throws: ZQEngineError.self) {
82+ _ = try await session.execute("margins")
83+ }
84+ await #expect(throws: ZQEngineError.self) {
85+ _ = try await session.execute("margins, dydx(nonexistent)")
86+ }
87+ }
88+}
modified Tests/Fixtures/generate.R +29 −0
@@ -291,6 +291,35 @@ meat_iv <- t(Xhat * as.vector(u_iv^2)) %*% Xhat * (Nw / df_iv)
291 V_iv_r <- XhXhinv %*% meat_iv %*% XhXhinv 291 V_iv_r <- XhXhinv %*% meat_iv %*% XhXhinv
292 emit("iv_se_hc1_price", sqrt(V_iv_r[1, 1])) 292 emit("iv_se_hc1_price", sqrt(V_iv_r[1, 1]))
293 293
294+# ---------------------------------------------------- margins (AME, delta)
295+# Average marginal effects with delta-method SEs at the converged
296+# coefficients, mirroring the engine's formulas exactly.
297+ame_glm <- function(fit, var, fprime, fsecond) {
298+ Xg <- model.matrix(fit)
299+ eta <- as.vector(Xg %*% coef(fit))
300+ fp <- fprime(eta)
301+ fs <- fsecond(eta)
302+ b <- coef(fit)[[var]]
303+ ame <- b * mean(fp)
304+ grad <- colMeans(Xg * fs) * b
305+ grad[[var]] <- grad[[var]] + mean(fp)
306+ V <- glm_bread(fit)
307+ se <- sqrt(as.numeric(t(grad) %*% V %*% grad))
308+ c(ame = ame, se = se)
309+}
310+
311+m_lg <- ame_glm(
312+ lg, "price",
313+ function(e) plogis(e) * (1 - plogis(e)),
314+ function(e) plogis(e) * (1 - plogis(e)) * (1 - 2 * plogis(e))
315+)
316+emit("margins_logit_price", m_lg[["ame"]])
317+emit("margins_logit_se", m_lg[["se"]])
318+
319+m_ps <- ame_glm(ps, "price", exp, exp)
320+emit("margins_pois_price", m_ps[["ame"]])
321+emit("margins_pois_se", m_ps[["se"]])
322+
294 # -------------------------------------------------------------- correlate 323 # -------------------------------------------------------------- correlate
295 emit("corr_rev_price", cor(d$revenue, d$price)) 324 emit("corr_rev_price", cor(d$revenue, d$price))
296 emit("corr_rev_logrev", cor(d$revenue, d$log_rev)) 325 emit("corr_rev_logrev", cor(d$revenue, d$log_rev))
297 326
Menu

Theme
↑↓navigate ↵open escclose /⌘Kanywhere