Add keyword arg to modelmatrix; define momentmatrix

src/regressionmodel.jl

@@ -34,19 +34,25 @@ Return the mean of the response.
 function meanresponse end
 
 """
-    modelmatrix(model::RegressionModel)
+    modelmatrix(model::RegressionModel; weighted::Bool=false)
 
-Return the model matrix (a.k.a. the design matrix).
+Return the model matrix (design matrix) or, if `weighted=true` the weighted 
+model matrix, i.e. `X * sqrt.(W)`, where `X` is the model matrix and
+`W` is the diagonal matrix whose elements are the [model weights](@ref weights(::StatisticalModel)). 
 """
-function modelmatrix end
+function modelmatrix(model::RegressionModel; weighted::Bool=false) end
 
 """
-    crossmodelmatrix(model::RegressionModel)
+    crossmodelmatrix(model::RegressionModel; weighted::Bool=false)
 
-Return `X'X` where `X` is the model matrix of `model`.
+Return `X'X` where `X` is the model/design matrix of `model` or, if `weighted=true`, `X'WX`, 
+where `W` is the diagonal matrix whose elements are the [model weights](@ref weights(::StatisticalModel)). 
 This function will return a pre-computed matrix stored in `model` if possible.
 """
-crossmodelmatrix(model::RegressionModel) = (x = modelmatrix(model); Symmetric(x' * x))
+function crossmodelmatrix(model::RegressionModel; weighted::Bool=false) 
+    x = weighted ? modelmatrix(model; weighted=weighted) : modelmatrix(model)
+    return Symmetric(x' * x)
+end
 
 """
     leverage(model::RegressionModel)
@@ -65,11 +71,13 @@ of each data point.
 function cooksdistance end
 
 """
-    residuals(model::RegressionModel)
+    residuals(model::RegressionModel; weighted::Bool=false)
+
+Return the residuals of the model or, if `weighted=true`, the residuals multiplied by
+the square root of the [model weights](@ref weights(::StatisticalModel)). 
 
-Return the residuals of the model.
 """
-function residuals end
+function residuals(model::RegressionModel; weighted::Bool=false) end
 
 """
     predict(model::RegressionModel, [newX])

src/statisticalmodel.jl

@@ -314,3 +314,24 @@ function adjr2(model::StatisticalModel, variant::Symbol)
 end
 
 const adjr² = adjr2
+
+"""
+    momentmatrix(model::StatisticalModel)
+
+Return the matrix containing the empirical moments defining the estimated parameters.
+
+For likelihood-based models, `momentmatrix` returns the log-score, i.e. the gradient
+of the log-likelihood evaluated at each observation. For semiparametric models, each row
+of the ``(n \\times k)`` matrix returned by `momentmatrix` is ``m(Z_i, b)``, where `m`
+is a vector-valued function evaluated at the estimated parameter `b` and ``Z_i`` is the 
+vector of data for entity `i`. The vector-valued function satisfies ``\\sum_{i=1}^n m(Z_i, b) = 0``.
+
+For linear and generalized linear models, the parameters of interest are the coefficients
+of the linear predictor. The moment matrix of a linear model is given by `u.*X`,
+where `u` is the vector of residuals and `X` is the model matrix. The moment matrix of
+a generalized linear model with link function `g`, the moment matrix is `X'e`, where `e`
+is given by ``Y-g^{-1}(X'b)`` where `X` is the model matrix, `Y` is the model response, and
+`b` is the vector of estimated coefficients. 
+"""
+function momentmatrix end
+

test/regressionmodel.jl

@@ -6,13 +6,31 @@ using StatsAPI: RegressionModel, crossmodelmatrix
 struct MyRegressionModel <: RegressionModel
 end
 
+struct MyWeightedRegressionModel <: RegressionModel
+    wts::AbstractVector
+end
+
 StatsAPI.modelmatrix(::MyRegressionModel) = [1 2; 3 4]
 
+function StatsAPI.modelmatrix(r::MyWeightedRegressionModel; weighted::Bool=false)
+    X = [1 2; 3 4]
+    weighted ? sqrt.(r.wts).*X : X
+end
+
+w = [0.3, 0.2]
+
 @testset "TestRegressionModel" begin
     m = MyRegressionModel()
+    r = MyWeightedRegressionModel(w)
 
     @test crossmodelmatrix(m) == [10 14; 14 20]
+    @test crossmodelmatrix(m; weighted=false) == [10 14; 14 20]
     @test crossmodelmatrix(m) isa Symmetric
+
+    @test crossmodelmatrix(r) == [10 14; 14 20]
+    @test crossmodelmatrix(r; weighted=false) == [10 14; 14 20]
+    @test crossmodelmatrix(r; weighted=true) ≈ [2.1 3.0; 3.0 4.4]
+    @test crossmodelmatrix(r; weighted=true) isa Symmetric
 end
 
 end # module TestRegressionModel

test/statisticalmodel.jl

@@ -36,4 +36,4 @@ StatsAPI.nobs(::MyStatisticalModel) = 100
     @test adjr2 === adjr²
 end
 
-end # module TestStatisticalModel
+end # module TestStatisticalModel