constrain<- {lava}

Add non-linear constraints to latent variable model

Description

Add non-linear constraints to latent variable model

Usage

## Default S3 replacement method:
constrain(x,par,args,endogenous=TRUE,...) <- value

## S3 replacement method for class '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

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.

Value

A lvm object.

Author(s)

Klaus K. Holst

See Also

regression, intercept, covariance

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)
## Likelihood of restricted model (MLE)
logLik(e1)

##############################
### 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))
constraints(e,data=d)
## 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)

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