spb/metrika Public

1//2//  GLMTests.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Foundation11import Testing12import ZQData13import ZQStats1415/// GLM estimators vs R glm() fixtures (epsilon 1e-12) at 1e-10 relative16/// tolerance.17@Suite("ZQGLM vs R fixtures")18struct GLMTests {19    let fixtures: Fixtures20    let purchase: [Double]21    let orders: [Double]22    let price: [Double]23    let firmID: [Int]2425    init() async throws {26        self.fixtures = try Fixtures()27        let store = try ZQDataStore()28        let frame = try await store.load(contentsOf: fixtures.datasetURL)29        let (purchaseAll, _) = try frame.requireNumeric("purchase")30        let (ordersAll, _) = try frame.requireNumeric("orders")31        let (priceAll, priceMissing) = try frame.requireNumeric("price")32        let (firmAll, _) = try frame.requireNumeric("firm_id")3334        var purchase: [Double] = [], orders: [Double] = []35        var price: [Double] = [], firm: [Int] = []36        for i in 0..<frame.rowCount where !priceMissing[i] {37            purchase.append(purchaseAll[i])38            orders.append(ordersAll[i])39            price.append(priceAll[i])40            firm.append(Int(firmAll[i]))41        }42        self.purchase = purchase43        self.orders = orders44        self.price = price45        self.firmID = firm46    }4748    @Test("logit matches R glm(binomial)")49    func logit() throws {50        let result = try ZQGLM.fit(51            y: purchase, predictors: [("price", price)], family: .logit52        )53        expectClose(result.coefficients[0].estimate, fixtures["logit_b_price"], "b[price]")54        expectClose(result.coefficients[1].estimate, fixtures["logit_b_cons"], "b[_cons]")55        expectClose(56            result.coefficients[0].standardError, fixtures["logit_se_price"], "se[price]"57        )58        expectClose(59            result.coefficients[1].standardError, fixtures["logit_se_cons"], "se[_cons]"60        )61        expectClose(result.logLikelihood, fixtures["logit_ll"], "log likelihood")62        expectClose(result.nullLogLikelihood, fixtures["logit_ll0"], "null log likelihood")63        expectClose(try #require(result.chiSquared), fixtures["logit_chi2"], "LR chi2")64        expectClose(65            try #require(result.chiSquaredPValue),66            fixtures["logit_chi2p"], rtol: 1e-9, "Prob > chi2"67        )68    }6970    @Test("logit robust and cluster standard errors")71    func logitRobust() throws {72        let robust = try ZQGLM.fit(73            y: purchase, predictors: [("price", price)], family: .logit, variance: .hc074        )75        expectClose(76            robust.coefficients[0].standardError,77            fixtures["logit_se_hc0_price"], "robust se[price]"78        )79        let clustered = try ZQGLM.fit(80            y: purchase, predictors: [("price", price)], family: .logit,81            variance: .cluster(firmID)82        )83        expectClose(84            clustered.coefficients[0].standardError,85            fixtures["logit_se_cluster_price"], "cluster se[price]"86        )87    }8889    @Test("probit matches R glm(binomial(probit))")90    func probit() throws {91        let result = try ZQGLM.fit(92            y: purchase, predictors: [("price", price)], family: .probit93        )94        expectClose(result.coefficients[0].estimate, fixtures["probit_b_price"], "b[price]")95        expectClose(result.coefficients[1].estimate, fixtures["probit_b_cons"], "b[_cons]")96        expectClose(97            result.coefficients[0].standardError, fixtures["probit_se_price"], "se[price]"98        )99        expectClose(result.logLikelihood, fixtures["probit_ll"], "log likelihood")100    }101102    @Test("poisson matches R glm(poisson)")103    func poisson() throws {104        let result = try ZQGLM.fit(105            y: orders, predictors: [("price", price)], family: .poisson106        )107        expectClose(result.coefficients[0].estimate, fixtures["pois_b_price"], "b[price]")108        expectClose(result.coefficients[1].estimate, fixtures["pois_b_cons"], "b[_cons]")109        expectClose(110            result.coefficients[0].standardError, fixtures["pois_se_price"], "se[price]"111        )112        expectClose(result.logLikelihood, fixtures["pois_ll"], "log likelihood")113114        let robust = try ZQGLM.fit(115            y: orders, predictors: [("price", price)], family: .poisson, variance: .hc0116        )117        expectClose(118            robust.coefficients[0].standardError,119            fixtures["pois_se_hc0_price"], "robust se[price]"120        )121    }122123    @Test("logit rejects non-binary outcomes")124    func binaryValidation() {125        #expect(throws: ZQStatsError.self) {126            _ = try ZQGLM.fit(127                y: orders, predictors: [("price", price)], family: .logit128            )129        }130    }131132    @Test("correlate matches R cor()")133    func correlate() async throws {134        let store = try ZQDataStore()135        let frame = try await store.load(contentsOf: fixtures.datasetURL)136        let (revenueAll, _) = try frame.requireNumeric("revenue")137        let (priceAll, priceMissing) = try frame.requireNumeric("price")138        var revenue: [Double] = [], priceComplete: [Double] = []139        for i in 0..<frame.rowCount where !priceMissing[i] {140            revenue.append(revenueAll[i])141            priceComplete.append(priceAll[i])142        }143        let matrix = ZQCorrelate.matrix(columns: [revenue, priceComplete])144        expectClose(matrix[0][1], fixtures["corr_rev_price"], "cor(revenue, price)")145        #expect(matrix[0][0] == 1 && matrix[1][1] == 1)146    }147148    @Test("chi-square CDF and normal quantile vs R")149    func chiSquare() {150        expectClose(151            ZQDistributions.chiSquareCDF(3.8, df: 1),152            fixtures["dist_pchisq_3p8_1"], rtol: 1e-12, "pchisq(3.8, 1)"153        )154        expectClose(155            ZQDistributions.chiSquareCDF(25, df: 4),156            fixtures["dist_pchisq_25_4"], rtol: 1e-12, "pchisq(25, 4)"157        )158        expectClose(159            ZQDistributions.normalQuantile(0.975),160            fixtures["dist_qnorm_0p975"], rtol: 1e-12, "qnorm(0.975)"161        )162    }163}164