leonardv / restriktor Goto Github PK

The parametric and model-based bootstrapped p-value for shifted/relocated closed convex cones is incorrect. The solution is to shift y and \theta. This will be fixed asap.

A small issue found from JASP, probably a corner case:

library(restriktor)

fit <- lm(mpg ~ vs + cyl, data = mtcars)
restriktor::goric(fit, "vs ~ cyl < 0")

> Error in lav_syntax_parse_rhs(rhs = rhs.formula[[2L]], op = op) : 
  lavaan ERROR: I'm confused parsing this line: <cyl0

Would be great it the error message returned the problematic line without alterations... Thanks!

Maybe I have misunderstood the 'mixing weight' functions, con_weights_boot.
Considering the simple 2D model with correlation rho:

library(restriktor)
wts <- function(rho) {
 con_weights_boot(VCOV=matrix(c(1,rho,rho,1),nrow=2,ncol=2),
   Amat=diag(2),meq=0L,R=9999L)
}
set.seed(1234)
wts(0)
wts(0.9995)

This is restriktor 0.1-80.911
restriktor is BETA software! Please report any bugs.
[1] 0.2474247 0.5031503 0.2494249
[1] 0.00590059 0.50085009 0.49324932

The first answer looks like what I would expect, namely choose(2,0:2) / 4, or
c(0.25,0.5,0.25).
However, for rho near 1, I expect the weights to go to c(0.5,0,0.5). That is
not observed in this case.

The package has been archived from CRAN due to the norm2 package being archived. tried to install an older version of norm2 from source but ran into fortran problems. Is there a way to work around this dependency?

https://cran.r-project.org/package=restriktor

Package ‘restriktor’ was removed from the CRAN repository.

Formerly available versions can be obtained from the [archive](https://cran.r-project.org/src/contrib/Archive/restriktor).

Archived on 2023-08-25 as requires archived package 'norm2'.

A summary of the most recent check results can be obtained from the [check results archive](https://cran-archive.r-project.org/web/checks/2023/2023-08-25_check_results_restriktor.html).

Please use the canonical form https://cran.r-project.org/package=restriktor to link to this page.

It is possible to obtain different results for model specifications that should be intuitively the same.

First, I tried to specify an equality constraint between three groups like this:

fit <- lm(Petal.Width~0+Species, iris)
con0 <- restriktor::restriktor(fit, constraints = "Speciessetosa == Speciesversicolor == Speciesvirginica")

however the analysis fails to return correct coefficients:

> coef(con0)
    Speciessetosa Speciesversicolor  Speciesvirginica 
            1.000             1.326             2.026

Which lead me to believe that the package does not handle multiple constraints within one statement properly.

The alternative how to specify this model would be a series of equality constraints:

con1 <- restriktor::restriktor(fit, constraints = "Speciessetosa == Speciesversicolor; Speciesversicolor == Speciesvirginica")

but we may also add another equality constraint (that is arguably redundant):

con2 <- restriktor::restriktor(fit, constraints = "Speciessetosa == Speciesversicolor; Speciesversicolor == Speciesvirginica; Speciessetosa == Speciesvirginica")

Although con1 and con2 seem to get correct coefficients, the model comparison suggests that con2 is better, because for some reason it assigns a smaller penalty to that model:

> restriktor::goric(con1, con2, comparison = "none")
restriktor (0.3-100): generalized order-restriced information criterion: 

Results:
   model    loglik  penalty    goric  goric.weights
1   con1  -171.615    2.000  347.229          0.269
2   con2  -171.615    1.000  345.229          0.731
---
Restriktor message: Since the constraint matrix is not full row-rank, the level probabilities 
 are calculated using mix.weights = "boot" (the default is mix.weights = "pmvnorm").
 For more information see ?restriktor.

This does not seem to be limited only to equality constraints, but also occurs in inequality constraints. Using the equality constraints is only easier for this example.

When calling the goric() function with type = 'goricc' and comparison = 'complement', the result-column with the goricc values is called goric rather than goricc. Although this is not a very big deal, it breaks my simulation in which I want to extract that specific column for both unconstrained and complement alternative hypotheses. I think that it can be solved by getting

names(df)[4] <- type

outside of the if (comparison == 'unconstrained') block on line 502, or to add it to all three if blocks

if (comparison == 'unconstrained') # line 502 in `con_goric.R`, already contains this code
if (comparison == 'none')          # line 550 in `con_goric.R`
if (comparison == 'complement')    # line 567 in `con_goric.R`

Good afternoon!

I encountered this one problem: I am estimating ARDL and ECM models, so the usual inference is done via sandwich::vcovHAC. Then, some linear hypotheses about the long-run elasticities are tested (equality constraints), and this information is incorporated into the model.

In theory, having an equality constraint is equivalent to rewriting the RHS of the equation with fewer parameters and subtracting the known information from both sides (like using contrasts). Then, given the variance-covariance of the reduced-rank model, one can obtain the standard errors for the original regressors via linear transformations. In principle, if one has a HAC estimator of the reduced-rank model variance, then, it can be used to go obtain the standard errors in the original model.

However, I have been unable to find a way to apply custom VCOV estimators because it is not trivial how to get the reduced model in the first place. It is obvious that the Kronecker product for yVx and Dmat are used, but then, we never see the model matrix with constraints, which is necessary to do the conversion.

Back in the day, I implemented a function (for over-identified GMM) that uses the reduced row-echelon form to create a mapping for the moment function given the constraint matrix. The user would specify the full model, supply the constraints, and the function would create a moment function incorporation the equality constraints. E.g. if there are 14 coefficients in a linear model and we have linear constraints for 2 of them (the matrix A is 2x14), then, we minimise the GMM criterion w.r.t. a parameter vector of length 12 that gets unwrapped into the proper RHS relationships using the reduced row-echelon form function mapping R^{12} |--> R^{14}. At the end of the day, we have 12 parameters and 14 orthogonality conditions, and we are happy with the GMM solution. Given the values of the 14-dimensional moment function in the optimum (where the FOC holds for 12 parameters), it is possible to construct any VCOV estimator (HC, HAC, cluster-robust etc.). In restriktor, it is trivial to get the residuals of the constrained model, but then, the sandwich matrix must not be constructed with the unconstrained model matrix (because the FOCs for the unconstrained OLS criterion minimisation do not hold at the constrained optimum). Instead, the meat and bread functions of the sandwich should operate on the model matrix of the constrained model, and the mapping back to the original variables should be linear.

Is there a way to apply vcovHAC to the constrained model and get autocorrelation-robust standard errors for it? Or is there any other way to get HAC-robust standard errors with the existing facilities of restriktor? If not, I may attempt adjusting this linear mapping for use with equality constraints.

Hi,

Thanks a lot for your work on this package.
I faced the could not find function "lqs" error while fitting a constraint RLM (using restriktor::restriktor(...)). It seems that MASS is not loaded in the environment prior to the instruction do.call("lqs", c(list(x, y, intercept = FALSE), lqs.control)) in conRLM_fit.R (line 79).
A straightforward workaround is to load MASS myself before fitting the model.

Have a nice day,
Vadim