## S3 method for class 'default':
constrain(x,par,args,...) <- value## S3 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,...)
lvm-objectargs argument).parlvm-objectlvm object.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.
############################## ### 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
[object Object]
regression, intercept,
covariance