###{{{ print.fix
##' @export
print.fix <- function(x,exo=FALSE,...) {
switch(attributes(x)$type,
reg = cat("Regression parameters:\n"),
cov = cat("Covariance parameters:\n"),
mean = cat("Intercept parameters:\n"))
M <- linconstrain(x,print=TRUE)
invisible(x)
}
linconstrain <- function(x,print=TRUE,indent=" ",exo=FALSE,...) {
idx <- seq_len(attributes(x)$nvar)
idx0 <- setdiff(idx,attributes(x)$exo.idx)
if (!exo & attributes(x)$type!="reg")
idx <- idx0
if (attributes(x)$type=="mean") {
if (length(idx)>0){
M <- rbind(unlist(x[idx]))
rownames(M) <- ""
M[is.na(M)] <- "*"
} else {
M <- NULL
}
} else {
if (length(x$rel)==0) {
M <- NULL
} else {
M <- x$rel[idx,idx,drop=FALSE]
M[M==0] <- NA
M[M==1] <- "*"
M[which(!is.na(x$labels[idx,idx]))] <- x$labels[idx,idx][which(!is.na(x$labels[idx,idx]))]
M[which(!is.na(x$values[idx,idx]))] <- x$values[idx,idx][which(!is.na(x$values[idx,idx]))]
if (attributes(x)$type=="reg")
M <- t(M[,idx0,drop=FALSE])
}
}
if (print) {
M0 <- M
if (NROW(M)>0)
rownames(M0) <- paste(indent,rownames(M))
print(M0,quote=FALSE,na.print="",...)
}
invisible(M)
}
###}}} print.fix
###{{{ intfix
##' @export
"intfix" <- function(object,...) UseMethod("intfix")
##' @export
"intfix<-" <- function(object,...,value) UseMethod("intfix<-")
##' Fix mean parameters in 'lvm'-object
##'
##' Define linear constraints on intercept parameters in a \code{lvm}-object.
##'
##'
##' The \code{intercept} function is used to specify linear constraints on the
##' intercept parameters of a latent variable model. As an example we look at
##' the multivariate regression model
##'
##' \deqn{ E(Y_1|X) = \alpha_1 + \beta_1 X} \deqn{ E(Y_2|X) = \alpha_2 + \beta_2
##' X}
##'
##' defined by the call
##'
##' \code{m <- lvm(c(y1,y2) ~ x)}
##'
##' To fix \eqn{\alpha_1=\alpha_2} we call
##'
##' \code{intercept(m) <- c(y1,y2) ~ f(mu)}
##'
##' Fixed parameters can be reset by fixing them to \code{NA}. For instance to
##' free the parameter restriction of \eqn{Y_1} and at the same time fixing
##' \eqn{\alpha_2=2}, we call
##'
##' \code{intercept(m, ~y1+y2) <- list(NA,2)}
##'
##' Calling \code{intercept} with no additional arguments will return the
##' current intercept restrictions of the \code{lvm}-object.
##'
##' @aliases intercept intercept<- intercept.lvm intercept<-.lvm intfix intfix
##' intfix<- intfix.lvm intfix<-.lvm
##' @param object \code{lvm}-object
##' @param vars character vector of variable names
##' @param value Vector (or list) of parameter values or labels (numeric or
##' character) or a formula defining the linear constraints (see also the
##' \code{regression} or \code{covariance} methods).
##' @param \dots Additional arguments
##' @usage
##' \method{intercept}{lvm}(object, vars, ...) <- value
##' @return
##'
##' A \code{lvm}-object
##' @note
##'
##' Variables will be added to the model if not already present.
##' @author Klaus K. Holst
##' @seealso \code{\link{covariance<-}}, \code{\link{regression<-}},
##' \code{\link{constrain<-}}, \code{\link{parameter<-}},
##' \code{\link{latent<-}}, \code{\link{cancel<-}}, \code{\link{kill<-}}
##' @keywords models regression
##' @export
##' @examples
##'
##'
##' ## A multivariate model
##' m <- lvm(c(y1,y2) ~ f(x1,beta)+x2)
##' regression(m) <- y3 ~ f(x1,beta)
##' intercept(m) <- y1 ~ f(mu)
##' intercept(m, ~y2+y3) <- list(2,"mu")
##' intercept(m) ## Examine intercepts of model (NA translates to free/unique paramete##r)
##'
##'
"intercept" <- function(object,...) UseMethod("intercept")
##' @export
intercept.lvm <- function(object,value,...) {
if (!missing(value)) {
intercept(object,...) <- value
return(object)
}
res <- object$mean; attr(res,"type") <- "mean"
attr(res,"exo.idx") <- index(object)$exo.idx
attr(res,"nvar") <- length(res)
class(res) <- "fix"
return(res)
}
##' @export
intfix.lvm <- intercept.lvm
##' @export
"intercept<-" <- function(object,...,value) UseMethod("intercept<-")
##' @export
"intercept<-.lvm" <- function(object, vars,...,value) {
if (!missing(vars) && inherits(value,"formula")) value <- all.vars(value)
if (inherits(value,"formula")) {
lhs <- getoutcome(value)
yy <- decomp.specials(lhs)
if ((inherits(value[[3]],"logical") && is.na(value[[3]]))) {
intfix(object,yy) <- NA
return(object)
}
tt <- terms(value)
xf <- attributes(terms(tt))$term.labels
res <- lapply(xf,decomp.specials)[[1]]
myvalue <- char2num(as.list(res))
myvalue <- lapply(myvalue, function(x) ifelse(x=="NA",NA,x))
intfix(object,yy) <- myvalue
object$parpos <- NULL
return(object)
}
if (inherits(vars,"formula")) {
vars <- all.vars(vars)
}
object$mean[vars] <- value
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
##' @export
"intfix<-.lvm" <- function(object, ..., value) {
intercept(object, ...) <- value
object
}
###}}} intfix
###{{{ covfix
##' @export
"covfix" <- function(object,...) UseMethod("covfix")
##' @export
covfix.lvm <- function(object,...) {
res <- list(rel=object$cov, labels=object$covpar, values=object$covfix); attr(res,"type") <- "cov"
attr(res,"exo.idx") <- index(object)$exo.idx
attr(res,"nvar") <- NROW(res$rel)
class(res) <- "fix"
return(res)
}
##' @export
"covfix<-" <- function(object,...,value) UseMethod("covfix<-")
##' @export
"covfix<-.lvm" <- function(object, var1, var2=var1, pairwise=FALSE, exo=FALSE, ..., value) {
if (inherits(var1,"formula")) {
var1 <- all.vars(var1)
}
if (inherits(var2,"formula")) {
var2 <- all.vars(var2)
}
object <- addvar(object,c(var1,var2),reindex=FALSE,...)
allvars <- c(var1,var2)
xorg <- exogenous(object)
exoset <- setdiff(xorg,allvars)
if (!exo & length(exoset)<length(xorg)) {
exogenous(object) <- exoset
}
if (inherits(value,"formula")) value <- all.vars(value)
if (pairwise) {
p <- 0
K <- length(var1)*(length(var1)-1)/2
if (length(value)==1)
value <- rep(value,K)
if (length(value)!=K) stop("Wrong number of parameters")
for (i in seq_len(length(var1)-1)) {
for (j in seq(i+1,length(var1))) {
p <- p+1
valp <- char2num(value[[p]])
if (is.na(value[[p]]) | value[[p]]=="NA") {
object$covfix[var1[i],var1[j]] <- object$covpar[var1[i],var1[j]] <- NA
object$covfix[var1[j],var1[i]] <- object$covpar[var1[j],var1[i]] <- NA
}
else {
object$cov[var1[i],var1[j]] <- object$cov[var1[j],var1[i]] <- 1
if (is.numeric(value[[p]]) | !is.na(valp)) {
object$covfix[var1[i],var1[j]] <- object$covfix[var1[j],var1[i]] <- valp
object$covpar[var1[i],var1[j]] <- object$covpar[var1[j],var1[i]] <- NA
} else {
object$covpar[var1[i],var1[j]] <- object$covpar[var1[j],var1[i]] <- value[[p]]
object$covfix[var1[i],var1[j]] <- object$covfix[var1[j],var1[i]] <- NA
}
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
if (is.null(var2)) {
if (length(value)==1)
value <- rep(value,length(var1))
if (length(value)!=length(var1)) stop("Wrong number of parameters")
for (i in seq_along(var1)) {
vali <- char2num(value[[i]])
if (is.na(value[[i]]) | value[[i]]=="NA") {
object$covfix[var1[i],var1[i]] <- object$covpar[var1[i],var1[i]] <- NA
}
else {
if (is.numeric(value[[i]]) | !is.na(vali)) {
object$covfix[var1[i],var1[i]] <- vali
object$covpar[var1[i],var1[i]] <- NA
} else {
object$covfix[var1[i],var1[i]] <- NA
object$covpar[var1[i],var1[i]] <- value[[i]]
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
if (length(var1)==length(var2) & length(var1)==length(value)) {
p <- 0
for (i in seq_along(var1)) {
p <- p+1
valp <- char2num(value[[p]])
if (is.na(value[[p]]) | value[[p]]=="NA") {
object$covfix[var1[i],var2[i]] <- object$covpar[var1[i],var2[i]] <- NA
object$covfix[var2[i],var1[i]] <- object$covpar[var2[i],var1[i]] <- NA
}
else {
object$cov[var1[i],var2[i]] <- object$cov[var2[i],var1[i]] <- 1
if (is.numeric(value[[p]]) | !is.na(valp)) {
object$covfix[var1[i],var2[i]] <- object$covfix[var2[i],var1[i]] <- valp
object$covpar[var1[i],var2[i]] <- object$covpar[var2[i],var1[i]] <- NA
} else {
object$covpar[var1[i],var2[i]] <- object$covpar[var2[i],var1[i]] <- value[[p]]
object$covfix[var1[i],var2[i]] <- object$covfix[var2[i],var1[i]] <- NA
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
K <- length(var1)*length(var2)
if (length(value)==1)
value <- rep(value,K)
if (length(value)!=K) stop("Wrong number of parameters")
p <- 0
for (i in seq_along(var1)) {
for (j in seq_along(var2)) {
if (!pairwise | var1[i]!=var2[j]) {
p <- p+1
valp <- char2num(value[[p]])
if (is.na(value[[p]]) | value[[p]]=="NA") {
object$covfix[var1[i],var2[j]] <- object$covpar[var1[i],var2[j]] <- NA
object$covfix[var2[j],var1[i]] <- object$covpar[var2[j],var1[i]] <- NA
}
else {
object$cov[var1[i],var2[j]] <- object$cov[var2[j],var1[i]] <- 1
if (is.numeric(value[[p]]) | !is.na(valp)) {
object$covfix[var1[i],var2[j]] <- object$covfix[var2[j],var1[i]] <- valp
object$covpar[var1[i],var2[j]] <- object$covpar[var2[j],var1[i]] <- NA
} else {
object$covpar[var1[i],var2[j]] <- object$covpar[var2[j],var1[i]] <- value[[p]]
object$covfix[var1[i],var2[j]] <- object$covfix[var2[j],var1[i]] <- NA
}
}
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
###}}} covfix
###{{{ regfix
##' @export
"regfix" <- function(object,...) UseMethod("regfix")
##' @export
regfix.lvm <- function(object,...) {
res <- list(rel=index(object)$M, labels=object$par, values=object$fix); attr(res,"type") <- "reg"
attr(res,"exo.idx") <- index(object)$exo.idx
attr(res,"nvar") <- NROW(res$rel)
class(res) <- "fix"
return(res)
}
##' @export
"regfix<-" <- function(object,...,value) UseMethod("regfix<-")
##' @export
"regfix<-.lvm" <- function(object, to, from=NULL, exo=lava.options()$exogenous, variance, y,x, ..., value) {
if (!missing(y)) {
if (inherits(y,"formula")) {
y <- all.vars(y)
}
to <- y
}
if (!missing(x)) {
if (inherits(x,"formula")) x <- all.vars(x)
from <- x
}
if (is.null(to)) stop("variable list needed")
if (inherits(to,"formula")) {
val <- procformula(object,to,exo=exo)
constrain_formula <- to
object <- val$object
ys <- val$ys
xs <- val$xs
if (!missing(variance))
covariance(object, ys) <- variance
to <- ys; from <- xs
} else {
object <- addvar(object, c(to, from), reindex=FALSE, ...)
newexo <- from
notexo <- to
curvar <- index(object)$var
if (exo) {
oldexo <- exogenous(object)
newexo <- setdiff(newexo, c(notexo, curvar))
exogenous(object) <- union(newexo, setdiff(oldexo, notexo))
}
}
if (inherits(value,"formula")) value <- all.vars(value)
nonlinear_function <- NULL
if (is.function(value)) {
vfun <- value
aname <- setdiff(names(formals(value)), "...")
if (is.null(to)) { # formula was apparently given as ~outcome
to <- from[1]
from <- NULL
}
pname <- parameter(object)
if (length(from)==0) {
constrain_formula <- toformula(to, aname)
from <- setdiff(aname, pname)
nonlinear_function <- function(...) {
as.vector(do.call(vfun,
structure(as.list(as.data.frame(...)), names=aname)))
}
} else {
nonlinear_function <- function(...) as.vector(vfun(...))
##nonlinear_function <- vfun
from <- setdiff(from, pname)
}
object <- addvar(object, c(to, from))
value <- 0
}
if (length(from)==length(to) && length(from)==length(value)
&& is.null(nonlinear_function)) {
for (i in seq_along(from)) {
if (object$M[from[i], to[i]]==0) { ## Not adjacent! ##!isAdjacent(Graph(object), from[i], to[i])) {
object <- regression(object, to=to[i], from=from[i])
}
vali <- char2num(value[[i]])
if (is.na(value[[i]]) | value[[i]]=="NA") {
object$fix[from[i], to[i]] <- object$par[from[i], to[i]] <- NA
}
else {
if (is.numeric(value[[i]]) | !is.na(vali)) {
object$fix[from[i], to[i]] <- vali
object$par[from[i], to[i]] <- NA
} else {
object$par[from[i], to[i]] <- value[[i]]
object$fix[from[i], to[i]] <- NA
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
return(object)
}
for (i in from) {
for (j in to) {
if (object$M[i,j]==0) { ##!isAdjacent(Graph(object), i, j)) {
object <- regression(object, to=j, from=i)
}
}
}
K <- length(from)*length(to)
if (length(value)==1)
value <- rep(value, K)
if (length(value)!=K) stop("Wrong number of parameters")
for (j in seq_along(to)) {
for (i in seq_along(from)) {
p <- (j-1)*length(from) + i
valp <- char2num(value[[p]])
if (is.na(value[[p]]) | value[[p]]=="NA")
object$fix[from[i], to[j]] <- object$par[from[i], to[j]] <- NA
else {
if (is.numeric(value[[p]]) | !is.na(valp)) {
object$fix[from[i], to[j]] <- valp
object$par[from[i], to[j]] <- NA
} else {
object$par[from[i], to[j]] <- value[[p]]
object$fix[from[i], to[j]] <- NA
}
}
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
index(object) <- reindex(object)
if (!is.null(nonlinear_function)) {
constrain(object, constrain_formula, endogenous=FALSE, ...) <- nonlinear_function
}
return(object)
}
###}}} regfix
###{{{ parfix
##' @export
"parfix<-" <- function(x,...,value) UseMethod("parfix<-")
##' @export
"parfix<-.lvm" <- function(x,idx,...,value) {
parfix(x,idx,value,...)
}
##' @export
"parfix" <- function(x,...) UseMethod("parfix")
## m <- lvm(c(y[m:v]~b*x))
## constrain(m,b~a) <- base::identity
##' @export
parfix.lvm <- function(x,idx,value,fix=FALSE,...) {
object <- Model(x)
if (fix)
object <- fixsome(object)
if (length(idx)!=length(value))
value <- rep(value,length.out=length(idx))
value <- as.list(value)
I <- index(object)
V <- with(I, matrices2(Model(object), seq_len(npar.mean+npar+npar.ex)))
V$A[I$M0!=1] <- 0; V$P[I$P0!=1] <- 0
v.fix <- which(V$v%in%idx) ## Intercepts
vval <- V$v[v.fix]
v.ord <- match(vval,idx)
Pval <- V$P[V$P%in%idx] ## Variance/covariance
P.fix <- which(matrix(V$P%in%idx,nrow=nrow(V$P)),arr.ind=TRUE)
P.ord <- match(Pval,idx)
Aval <- V$A[which(V$A%in%idx)] ## Regression parameters
A.fix <- which(matrix(V$A%in%idx,nrow=nrow(V$A)),arr.ind=TRUE)
A.ord <- match(Aval,idx)
e.fix <- which(V$e%in%idx)
eval <- V$e[e.fix]
e.ord <- match(eval,idx)
for (i in seq_len(length(e.fix))) {
object$exfix[[e.fix[i]]] <- value[[e.ord[i]]]
}
for (i in seq_len(length(v.fix))) {
object$mean[[v.fix[i]]] <- value[[v.ord[i]]]
}
for (i in seq_len(nrow(A.fix))) {
if (is.numeric(value[[ A.ord[i] ]])){
object$fix[A.fix[i,1],A.fix[i,2]] <- value[[A.ord[i]]]
object$par[A.fix[i,1],A.fix[i,2]] <- NA
} else {
object$par[A.fix[i,1],A.fix[i,2]] <- value[[A.ord[i]]]
object$fix[A.fix[i,1],A.fix[i,2]] <- NA
}
}
for (i in seq_len(nrow(P.fix))) {
if (is.numeric(value[[ P.ord[i] ]])) {
object$covfix[P.fix[i,1],P.fix[i,2]] <- value[[P.ord[i]]]
object$covpar[P.fix[i,1],P.fix[i,2]] <- NA
} else {
object$covpar[P.fix[i,1],P.fix[i,2]] <- value[[P.ord[i]]]
object$covfix[P.fix[i,1],P.fix[i,2]] <- NA
}
}
newindex <- reindex(object)
object$parpos <- NULL
index(object)[names(newindex)] <- newindex
attributes(object)$fixed <- list(v=v.fix,A=A.fix,P=P.fix,e=e.fix)
return(object)
}
###}}} parfix
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.