Add non-linear constraints to latent variable model

# S3 method for default
constrain(x,par,args,endogenous=TRUE,...) <- value

# S3 method for multigroup
constrain(x,par,k=1,...) <- value

constraints(object,data=model.frame(object),vcov=object$vcov,level=0.95,
                        p=pars.default(object),k,idx,...)

Arguments

x

lvm-object

...

Additional arguments to be passed to the low level functions

value

Real function taking args as a vector argument

par

Name of new parameter. Alternatively a formula with lhs specifying the new parameter and the rhs defining the names of the parameters or variable names defining the new parameter (overruling the args argument).

args

Vector of variables names or parameter names that are used in defining par

endogenous

TRUE if variable is endogenous (sink node)

k

For multigroup models this argument specifies which group to add/extract the constraint

object

lvm-object

data

Data-row from which possible non-linear constraints should be calculated

vcov

Variance matrix of parameter estimates

level

Level of confidence limits

p

Parameter vector

idx

Index indicating which constraints to extract

Value

A lvm object.

Details

Add non-linear parameter constraints as well as non-linear associations between covariates and latent or observed variables in the model (non-linear regression).

As an example we will specify the follow multiple regression model:

$$E(Y|X_1,X_2) = \alpha + \beta_1 X_1 + \beta_2 X_2$$ $$V(Y|X_1,X_2) = v$$

which is defined (with the appropiate parameter labels) as

m <- lvm(y ~ f(x,beta1) + f(x,beta2))

intercept(m) <- y ~ f(alpha)

covariance(m) <- y ~ f(v)

The somewhat strained parameter constraint $$ v = \frac{(beta1-beta2)^2}{alpha}$$

can then specified as

constrain(m,v ~ beta1 + beta2 + alpha) <- function(x) (x[1]-x[2])^2/x[3]

A subset of the arguments args can be covariates in the model, allowing the specification of non-linear regression models. As an example the non-linear regression model $$ E(Y\mid X) = \nu + \Phi(\alpha + \beta X)$$ where \(\Phi\) denotes the standard normal cumulative distribution function, can be defined as

m <- lvm(y ~ f(x,0)) # No linear effect of x

Next we add three new parameters using the parameter assigment function:

parameter(m) <- ~nu+alpha+beta

The intercept of \(Y\) is defined as mu

intercept(m) <- y ~ f(mu)

And finally the newly added intercept parameter mu is defined as the appropiate non-linear function of \(\alpha\), \(\nu\) and \(\beta\):

constrain(m, mu ~ x + alpha + nu) <- function(x) pnorm(x[1]*x[2])+x[3]

The constraints function can be used to show the estimated non-linear parameter constraints of an estimated model object (lvmfit or multigroupfit). Calling constrain with no additional arguments beyound x will return a list of the functions and parameter names defining the non-linear restrictions.

The gradient function can optionally be added as an attribute grad to the return value of the function defined by value. In this case the analytical derivatives will be calculated via the chain rule when evaluating the corresponding score function of the log-likelihood. If the gradient attribute is omitted the chain rule will be applied on a numeric approximation of the gradient.

Author

Klaus K. Holst

Examples

##############################
### Non-linear parameter constraints 1
##############################
m <- lvm(y ~ f(x1,gamma)+f(x2,beta))
covariance(m) <- y ~ f(v)
d <- sim(m,100)
m1 <- m; constrain(m1,beta ~ v) <- function(x) x^2
## Define slope of x2 to be the square of the residual variance of y
## Estimate both restricted and unrestricted model
e <- estimate(m,d,control=list(method="NR"))
e1 <- estimate(m1,d)
p1 <- coef(e1)
p1 <- c(p1[1:2],p1[3]^2,p1[3])
## Likelihood of unrestricted model evaluated in MLE of restricted model
logLik(e,p1)
#> 'log Lik.' -144.784 (df=4)
## Likelihood of restricted model (MLE)
logLik(e1)
#> 'log Lik.' -144.784 (df=3)

##############################
### Non-linear regression
##############################

## Simulate data
m <- lvm(c(y1,y2)~f(x,0)+f(eta,1))
latent(m) <- ~eta
covariance(m,~y1+y2) <- "v"
intercept(m,~y1+y2) <- "mu"
covariance(m,~eta) <- "zeta"
intercept(m,~eta) <- 0
set.seed(1)
d <- sim(m,100,p=c(v=0.01,zeta=0.01))[,manifest(m)]
d <- transform(d,
               y1=y1+2*pnorm(2*x),
               y2=y2+2*pnorm(2*x))

## Specify model and estimate parameters
constrain(m, mu ~ x + alpha + nu + gamma) <- function(x) x[4]*pnorm(x[3]+x[1]*x[2])
 ## Reduce Ex.Timings
e <- estimate(m,d,control=list(trace=1,constrain=TRUE))
#>   0:     95.152055: -1.59936 -0.105361  0.00000  0.00000  0.00000
#>   1:     38.372165: -2.18276 0.313970  0.00000  0.00000 0.695563
#>   2:    -48.631398: -5.42107 -0.457717  1.00498 0.959271  2.42405
#>   3:    -105.61594: -4.91886 -0.355240  1.37309 0.129419  1.46668
#>   4:    -152.46013: -4.25655 -1.25021  1.67211 -0.125548  2.24731
#>   5:    -180.78974: -4.81503 -2.47787  1.62885 -0.253601  1.83976
#>   6:    -211.28334: -4.80871 -2.49722  1.70839 -0.00201505  2.15835
#>   7:    -221.43760: -4.99751 -2.78668  1.85800 -0.174273  2.15106
#>   8:    -234.66132: -4.89742 -3.07152  1.81950 0.0648691  2.00370
#>   9:    -235.81854: -4.84621 -3.39006  1.92584 -0.172017  2.00849
#>  10:    -243.54462: -4.84315 -3.41209  1.93840 -0.0450503  2.11931
#>  11:    -249.92699: -4.91487 -3.54656  1.94515 -0.0534665  2.04374
#>  12:    -255.87475: -4.87632 -3.87588  1.88862 -0.0315561  2.09476
#>  13:    -257.75915: -4.87427 -3.99772  1.89083 -0.0859948  2.04418
#>  14:    -260.62803: -4.95794 -4.02217  1.97157 -0.00862611  2.02770
#>  15:    -261.37043: -4.93303 -4.07131  1.96533 0.00792647  1.99765
#>  16:    -262.57217: -4.90751 -4.12615  1.96599 0.0129969  2.02148
#>  17:    -263.47532: -4.89103 -4.18616  1.97601 0.000431286  2.01046
#>  18:    -264.87912: -4.85970 -4.31134  1.98013 0.00245754  2.02871
#>  19:    -265.31715: -4.81146 -4.38793  1.88695 -0.00882105  2.02646
#>  20:    -266.24010: -4.85667 -4.48840  1.95621 -0.00917811  2.01806
#>  21:    -266.36904: -4.85545 -4.49006  1.95711 -0.00250353  2.02374
#>  22:    -266.42668: -4.85202 -4.49488  1.95787 0.000747066  2.01776
#>  23:    -266.52881: -4.83733 -4.50258  1.96390 0.00407995  2.02003
#>  24:    -266.85274: -4.78248 -4.55687  1.98001 0.00102902  2.01598
#>  25:    -267.11141: -4.79064 -4.63495  1.97188 0.00301230  2.01884
#>  26:    -267.25836: -4.76836 -4.70965  1.98157 0.00925814  2.01327
#>  27:    -267.27808: -4.74963 -4.75607  1.97690 0.00626678  2.01510
#>  28:    -267.27984: -4.74889 -4.75331  1.98103 0.00876767  2.01359
#>  29:    -267.28003: -4.75020 -4.75214  1.98027 0.00821352  2.01379
#>  30:    -267.28003: -4.75002 -4.75222  1.98016 0.00820164  2.01385
#>  31:    -267.28003: -4.75004 -4.75224  1.98020 0.00821070  2.01384
#>  32:    -267.28003: -4.75004 -4.75223  1.98020 0.00821063  2.01384
constraints(e,data=d)
#>    Estimate Std. Error  Z value Pr(>|z|)     2.5%    97.5%
#> mu 1.598142 0.02227334 71.75133        0 1.554487 1.641797
## Plot model-fit
plot(y1~x,d,pch=16); points(y2~x,d,pch=16,col="gray")
x0 <- seq(-4,4,length.out=100)
lines(x0,coef(e)["nu"] + coef(e)["gamma"]*pnorm(coef(e)["alpha"]*x0))



##############################
### Multigroup model
##############################
### Define two models
m1 <- lvm(y ~ f(x,beta)+f(z,beta2))
m2 <- lvm(y ~ f(x,psi) + z)
### And simulate data from them
d1 <- sim(m1,500)
d2 <- sim(m2,500)
### Add 'non'-linear parameter constraint
constrain(m2,psi ~ beta2) <- function(x) x
## Add parameter beta2 to model 2, now beta2 exists in both models
parameter(m2) <- ~ beta2
ee <- estimate(list(m1,m2),list(d1,d2),control=list(method="NR"))
summary(ee)
#> ||score||^2= 2.64252e-17 
#> Latent variables: 
#> ____________________________________________________
#> Group 1 (n=500)
#>                     Estimate Std. Error  Z value Pr(>|z|)
#> Regressions:                                             
#>    y~x               0.93621    0.04767 19.64130   <1e-12
#>    y~z               1.03489    0.03217 32.17115   <1e-12
#> Intercepts:                                              
#>    y                -0.05578    0.04814 -1.15869   0.2466
#> Residual Variances:                                      
#>    y                 1.15827    0.07326 15.81139         
#> ____________________________________________________
#> Group 2 (n=500)
#>                        Estimate Std. Error  Z value Pr(>|z|)
#> Regressions:                                                
#>    y~z                  0.96338    0.04278 22.51997   <1e-12
#> Intercepts:                                                 
#>    y                   -0.01359    0.04657 -0.29191   0.7704
#> Additional Parameters:                                      
#>    beta2                1.03489    0.03217 32.17115   <1e-12
#> Residual Variances:                                         
#>    y                    1.07306    0.06787 15.81139         
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Non-linear constraints:
#>     Estimate Std. Error  Z value      Pr(>|z|)      2.5%    97.5%
#> psi 1.034886 0.03216812 32.17115 4.470275e-227 0.9718373 1.097934
#> ────────────────────────────────────────────────────────────────────────────────
#> Estimator: gaussian 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#>  Number of observations = 1000 
#>  BIC = 2994.955 
#>  AIC = 2960.6 
#>  log-Likelihood of model = -1473.3 
#> 
#>  log-Likelihood of saturated model = -1473.128 
#>  Chi-squared statistic: q = 0.3442204 , df = 1 
#>   P(Q>q) = 0.5574032 
#> 
#>  RMSEA (90% CI): 0 (0;0.0698)
#>   P(RMSEA<0.05)=0.85508
#> 
#> rank(Information) = 7 (p=7)
#> condition(Information) = 5.215445
#> mean(score^2) = 3.775029e-18 
#> ────────────────────────────────────────────────────────────────────────────────

m3 <- lvm(y ~ f(x,beta)+f(z,beta2))
m4 <- lvm(y ~ f(x,beta2) + z)
e2 <- estimate(list(m3,m4),list(d1,d2),control=list(method="NR"))
e2
#> ____________________________________________________
#> Group 1 (n=500)
#>                     Estimate Std. Error  Z value Pr(>|z|)
#> Regressions:                                             
#>    y~x               0.93621    0.04767 19.64130   <1e-12
#>    y~z               1.03489    0.03217 32.17115   <1e-12
#> Intercepts:                                              
#>    y                -0.05578    0.04814 -1.15869   0.2466
#> Residual Variances:                                      
#>    y                 1.15827    0.07326 15.81139         
#> ____________________________________________________
#> Group 2 (n=500)
#>                     Estimate Std. Error  Z value Pr(>|z|)
#> Regressions:                                             
#>    y~x               1.03489    0.03217 32.17115   <1e-12
#>    y~z               0.96338    0.04278 22.51997   <1e-12
#> Intercepts:                                              
#>    y                -0.01359    0.04657 -0.29191   0.7704
#> Residual Variances:                                      
#>    y                 1.07306    0.06787 15.81139         
#>