Add keyword arg to modelmatrix; define momentmatrix

[gragusa]

1. The modelmatrix has now a keyword weighted=false which is useful for dealing with weighted models.
2. Add momentmatrix - this function is intended to return the matrix of estimating equations; for instance, for a linear model should return u*X, where u is the vector of residuals and X is the model matrix.

[nalimilan]

Suggested change

	wts
	wts::AbstractVector

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 75.00000% with 1 line in your changes missing coverage. Please review.

Project coverage is 97.43%. Comparing base (20b38e1) to head (93f8742).
Report is 9 commits behind head on main.

Files with missing lines	Patch %	Lines
src/regressionmodel.jl	75.00%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##              main      #16      +/-   ##
===========================================
- Coverage   100.00%   97.43%   -2.57%     
===========================================
  Files            3        2       -1     
  Lines           37       39       +2     
===========================================
+ Hits            37       38       +1     
- Misses           0        1       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[nalimilan]

Why transpose X? It sounds weird to change the orientation of the result depending on whether it's weighted or not.

[gragusa]

My bad...I will fix it

[mschauer]

Suggested change

	Return the model matrix (a.k.a. the design matrix) or, if `weighted=true` the weighted
	Return the model matrix (design matrix) or, if `weighted=true` the weighted

[mschauer]

Can we define weights?

[nalimilan]

We could add a link to the weights(::StatisticalModel) method. Indeed there can be confusion between prior weights and working weights (though these terms can also confuse casual users).

[gragusa]

How exactly do I add a link to weights(::StatisticalModel)? Is there a way to link docs from different packages?

[nalimilan]

It's in the same package so I think something like [model weights](@ref weights(::StatisticalModel)) should work. Better test it though by building the StatsBase docs (julia docs/make.jl) using the updated StatsAPI.

[gragusa]

The only thing that we probably should do is to allow for modelmatrix to take another keyword argument, e.g., droppcollinearcols, to return only the columns corresponding to non-NaN coefficients in a Pivot Cholesky.

We could do it next (after I drop a bomb-PR against GLM. The GLM PR is waiting for this PR to get merged)

[nalimilan]

Let's tackle this separately. :-)

I'd rather review the GLM PR before merging this one, usually having the implementation is a good way to check that the API is the right one.

[gragusa]

I think it would be helpful to think about the API for dealing with rank-deficient models. For instance,
‘modelmatrix’ returns all the columns even those related to coefficients that cannot be estimated. I think this is fine as a default, but it would be useful to have a keyword argument to return only the columns corresponding to the estimable coefficients and a mechanism to identify the indexes of these columns. I have hacked these second part in this PR.

[nalimilan]

Where does the square root come from exactly? Doesn't that assume a particular definition of residuals (i.e. using L2-norm rather than e.g. L1-norm)?

[gragusa]

Well, this is tricky. Also modelmatrix multiplies the entries of X by the square-root of the weights. Why?

Think about the linear model. With weights, the crossmodel matrix is $X'WX$. Then, to obtain it we can now do modelmatrix(lm1; weighted = true)'modelmatrix(lm1; weighted = true).

Notice that this is consistent with R; see, e.g., the function weighted.residuals which is in stats.

With weights, any weights is
$$\sqrt{w_i}y_i = \sqrt{w_i}x_i \beta + \sqrt{w_i}u_i.$$ So the understanding is that weighting single constituents of the model (y,x,u) amount to weight by $\sqrt{w_i}$.

[nalimilan]

Yeah, another tricky point. My understanding is that for residuals, the square root comes from the fact that deviance residuals themselves are defined as the square root of quantities which are partitions of the deviance. Right?

Note the R docstring for weighted.residuals says "Weighted residuals are based on the deviance residuals", which are only one kind of residual. Actually in R residuals also returns weighted residuals, except for response residuals, which are always unweighted. Maybe to be completely accurate we could say "for deviance and Pearson residuals...", so that packages are free to use different definitions (or throw an error) if needed?

[bkamins]

I think what @nalimilan says is that the assumption here (and in your change of modelmatrix) is that for all kinds of weights the weighted model matrix is X * sqrt.(W). Is it always true for FrequencyWeights, AnalyticWeights and ProbabilityWeights? x-ref: JuliaStats/GLM.jl#487

[gragusa]

@bkamins weighted residuals, weighted model matrix do not exist in statistics. They are only useful from a coding point of view - they make it easier to write neater code.

I have always defined these quantities as multiplied by $\sqrt{w_i}$ as it is much more convenient. Some thing for R — which returns silently squared-root weighted residuals. Also other packages, notable FixedEffectModels.jl does that.

@nalimilan make sense what you propose - I will add more context to the doc

[nalimilan]

Even if these don't exist in statistics, the question can be phrased as "are there situations where the returned value is useful, even when you don't know the kind of weights used". I think the answer is yes, but it's tricky, so... R base only supports analytic weights so it's not a great reference.

[gragusa]

Something that would be very useful and that I would like to add o this PR is invloglikhessian - inverse of the hessian of the likelihood. The name is not great, but invloglikelihoodhessian is a monster, and I don't like invloglikhess (but I would do anything to have this merged.

With this method defined (whose implementation for GLM is part of JuliaStats/GLM.jl#487) CovarianceMatrices.jl could drop the GLM dependency, and implementation of all the covariances could be done via the API.

[nalimilan]

Sure. Would it make sense to call it invhessian? Do other modeling packages define this function currently?

[gragusa]

Sure. Would it make sense to call it invhessian? Do other modeling packages define this function currently?

I don't think so. The R package dealing with robust variances uses meat and bread which is not very clear (but in that case does not matter since R is not composable and methods can be extended only by the package itself).

We already have the inverse of the log hessian of the likelihood in glm -- is invchol. For linear models, the inverse of the normal log-likelihood is (X'X)^{-1}. But for general likelihoods this is not the case.

Now, invhessian is a little too generic, but I could live with it.

[nalimilan]

I don't have a strong preference, but at least for consistency I think we should spell "likelihood" in full if we use that term. Luckily autocompletion will almost always work. ;-)