constrain<-: Add non-linear constraints to latent variable model

Description

Add non-linear constraints to latent variable model

Usage

# 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,...)

Value

A lvm object.

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

Author

Klaus K. Holst

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}) = α + β_{1} X_{1} + β_{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{(b e t a 1 - b e t a 2)^{2}}{a l p h a}$

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 ∣ X) = ν + Φ (α + β X)$ where $Φ$ 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 $α$ , $ν$ and $β$ :

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.

Examples

Run this code

##############################
### 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

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning