Ridge Option for non-positive definite sample covariance matrix

[Julia Walther]

Hi,

I have some questions and notes concerning the ridge option(s) that ought to be evoked when the sample covariance matrix is non-positive definite (in single-level SEM). Was not sure whether to post here or on github, since I am not sure whether I am using it wrong or it‘s a coding issue. Decided to try here first. Below are my questions/notes and some code.

Thanks in advance for your help!

1) What is the difference between the parameters "ridge" and "ridge.x"? How to use them?

Both are parameters of lavaan(), whose defaults we get with lavOptions(), but they are not explained in the doc*, see ?lavOptions.

I assume that while "ridge.x" ridges only exogenous variances, "ridge" ridges all variances. However, it seems that only ridge.x works. Though, I am not sure whether it works since fit@SampleStats@ridge, which I assumed to give a non-zero value when it was used, always gives 0. Or is it just for "ridge" and not "ridge.x", but then, why is there no fit@SampleStats@ridge.x? Or what does fit@SampleStats@ridge indicate? Is there a way to check whether "ridge" or "ridge.x" was actually applied?

*Actually "ridge" is explained in the doc but seems outdated. "ridge: Numeric. Small constant used for ridging. Only used if the sample covariance matrix is non positive definite.", but actually "ridge" is boolean and "ridge.constant" is numeric.

Two observations speak against the notion (or the proper performance) of "ridge" ridging all variances (see code below):

- In fit_1 to fit_4, where an undirected analysis model is used.

- In fit_5, where we have an directed model, ridge is also not evoked.

In general, why is ridge applied to the exogenous variances? To my knowledge, it‘s usually applied to all variances (e.g. QUELLEN).

2) What does the parameter "sample.icov" do?

It is not explained in ?lavOptions but looking at lav_samplestats_icov(), apparently the inverse of the sample covmat is computed by choleski decomposition. Since this requires a pd covmat, (exogenous) ridging is applied in case the covmat is npd.

However, why does fit_6@SampleStats@icov then does not give the same as lavaan:::inv.chol(S_pd.x) (i.e. trying the function manually with the exogenous ridged covmat). What does fit_6@SampleStats@icov yield?

Besides lav_samplestats_icov(), there is also lav_matrix_symmetric_inverse() that enables to find the inverse of a matrix (with more options; not only choleski decomp). When is why used? Why does fit_8 (i.e. when a pd covmat is given) only throw an error related to lav_matrix_symmetric_inverse() but not to lav_samplestats_icov()? Why can we turn "sample.icov" FALSE? Is then (only) the inverse found by lav_matrix_symmetric_inverse() used?

# in general: if lavaan ridge is evoked depends on model specification

#### 1) no ridge evoked when no exogenous variables (i.e. undirected relations only)

popModel <- "y ~~ 1* y; \n

x ~~ 1* x"

N <- 3

set.seed(123)

data <- simulateData(popModel, sample.nobs=N, model.type="lavaan")

rcx <- 1e-5 # lavaan ridge value, see lavOptions()$ridge.constant.x, still used to ridge all variances in the following

# create npd S by substracting lavaan ridge value from all variances:

S <- cov(data)

S_npd <- S

while (is.positive.definite(S_npd)) {

diag(S_npd) <- diag(S_npd) - rcx

}

# test if ridging all variances once (i.e., assumed ridge behavior) will make S_npd pd:

is.positive.definite(S_npd) # F

S_pd <- S_npd

diag(S_pd) <- diag(S_pd) + rcx

is.positive.definite(S_pd) # indeed T again

corrModel <- "y ~~ y; \n # null-modell

x ~~ x; \n

y ~~ x"

# only error received when S npd, ridge or ridge.x not evoked apparently:

fit_1 <- lavaan(corrModel,

sample.cov = S_npd,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge=T,

ridge.constant=rcx)

fit_2 <- lavaan(corrModel,

sample.cov = S_npd,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge.x=T)

# but when S is pd, ridge and ridge.x can be turned on however:

fit_3 <- lavaan(corrModel,

sample.cov = S,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge=T, # works, (i.e., fit@Options[["ridge"]] set to T) but not necessary here

ridge.constant=rcx)

fit_4 <- lavaan(corrModel,

sample.cov = S,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge.x=T) # works, (i.e., fit@Options[["ridge.x"]] set to T) but not necessary here

### 2) ridge evoked when regression model specified (i.e. directed relations present) (?)

regModel <- "y ~ x"

# fit object created but not sure whether ridge was actually applied

fit_5 <- lavaan(regModel,

sample.cov = S_npd,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge=T,

ridge.constant=rcx)

fit_5@SampleStats@ridge # 0, meaning not used??

# create S that is npd by substracting lavaan ridge value from exogenous variable:

S_npd.x <- S

while (is.positive.definite(S_npd.x)) {

diag(S_npd.x)[2] <- diag(S_npd.x)[2] - rcx

}

# test if ridging exogenous variables once (i.e. ridge.x behavior) will make S_npd.x pd:

is.positive.definite(S_npd.x) # F

S_pd.x <- S_npd.x

diag(S_pd.x) <- diag(S_pd.x) + rcx

is.positive.definite(S_pd.x) # indeed T again

fit_6 <- lavaan(regModel,

sample.cov = S_npd.x,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge=T,

ridge.constant=rcx)

# still error thrown, as before

fit_6@SampleStats@cov[[1]] == S_npd.x # input (i.e. not ridged)

fit_6@SampleStats@ridge # 0, meaning not used?

fit_6@Fit@converged # F

fit_6@SampleStats@x.idx # index of exogenous variables

fit_6@SampleStats@cov.x # variance of exogenous variable (i.e., what would get ridged), equals S_npd.x (i.e., not ridged)

fit_6@SampleStats@icov # inverted S?

lavaan:::inv.chol(S_npd.x) # does not work bc npd

lavaan:::inv.chol(S_pd.x) # works, but differs from fit_6@SampleStats@icov

# what changes when sample.icov=F?

fit_7 <- lavaan(regModel,

sample.cov = S_npd.x,

sample.nobs = N,

sample.cov.rescale=FALSE,

sample.icov=F, #

ridge=T,

ridge.constant=rcx)

fit_6@Options[["sample.icov"]] # T

fit_7@Options[["sample.icov"]] # F

# likelihood still same:

fit_6@loglik[["loglik"]]

fit_7@loglik[["loglik"]]

# ought to work since covmat is pd, but does not (inverse not found)

fit_8 <- lavaan(regModel,

sample.cov = S_pd.x,

sample.nobs = N,

sample.cov.rescale=FALSE,

ridge=T,

ridge.constant=rcx)

# View(lavaan:::lav_samplestats_icov)

# View(lavaan:::lav_matrix_symmetric_inverse)

[Julia Walther]

Well, forgot to cite, so here it is:

Yuan, K.-H., Wu, R., & Bentler, P. M. (2011). Ridge structural equation modelling with correlation matrices for ordinal and continuous data: Ridge SEM with correlation matrices. British Journal of Mathematical and Statistical Psychology, 64(1), 107–133. https://doi.org/10/cwd74t

[Yves Rosseel]

Hello Julia,

> 1) What is the difference between the parameters "ridge" and "ridge.x"?
> How to use them?

The documentation AND code for 'ridge' needs updating/fixing indeed. But
the ridge.x option is not used, and should have been removed.

Ridging (if requested) is currently only done for exogenous variables.
The goal was (only) to allow for rank-deficient design matrices. It was
never designed for ridge SEM (as in Yuan & Bentler, 2011).

However, as I just found out, it does not work if you provide the sample
covariance matrix directly! This is a bug that will be fixed soon.

> 2) What does the parameter "sample.icov" do?

If set to FALSE, it avoids the computation of the inverse of the sample
covariance matrix. This is only used internally (in sam()) to handle
special cases where the number of variables is smaller than the number
of observations. It is not supposed to be a user-visible option.

Yves.

[Yves Rosseel]

The ridge= option has been restored in the github version of lavaan.
Using ridge = TRUE will (as documented) add a small constant to all the
elements of the diagonal of the sample covariance matrix.

Yves.

[Julia Walther]

Thank you so much for your help, Yves!

Best,
Julia