1
2##' Add regression association to latent variable model
3##'
4##' Define regression association between variables in a \code{lvm}-object and
5##' define linear constraints between model equations.
6##'
7##'
8##' The \code{regression} function is used to specify linear associations
9##' between variables of a latent variable model, and offers formula syntax
10##' resembling the model specification of e.g. \code{lm}.
11##'
12##' For instance, to add the following linear regression model, to the
13##' \code{lvm}-object, \code{m}:
14##' \deqn{ E(Y|X_1,X_2) = \beta_1 X_1 + \beta_2 X_2}
15##' We can write
16##'
17##' \code{regression(m) <- y ~ x1 + x2}
18##'
19##' Multivariate models can be specified by successive calls with
20##' \code{regression}, but multivariate formulas are also supported, e.g.
21##'
22##' \code{regression(m) <- c(y1,y2) ~ x1 + x2}
23##'
24##' defines
25##' \deqn{ E(Y_i|X_1,X_2) = \beta_{1i} X_1 + \beta_{2i} X_2 }
26##'
27##' The special function, \code{f}, can be used in the model specification to
28##' specify linear constraints. E.g. to fix \eqn{\beta_1=\beta_2}
29##' , we could write
30##'
31##' \code{regression(m) <- y ~ f(x1,beta) + f(x2,beta)}
32##'
33##' The second argument of \code{f} can also be a number (e.g. defining an
34##' offset) or be set to \code{NA} in order to clear any previously defined
35##' linear constraints.
36##'
37##' Alternatively, a more straight forward notation can be used:
38##'
39##' \code{regression(m) <- y ~ beta*x1 + beta*x2}
40##'
41##' All the parameter values of the linear constraints can be given as the right
42##' handside expression of the assigment function \code{regression<-} (or
43##' \code{regfix<-}) if the first (and possibly second) argument is defined as
44##' well. E.g:
45##'
46##' \code{regression(m,y1~x1+x2) <- list("a1","b1")}
47##'
48##' defines \eqn{E(Y_1|X_1,X_2) = a1 X_1 + b1 X_2}. The rhs argument can be a
49##' mixture of character and numeric values (and NA's to remove constraints).
50##'
51##' The function \code{regression} (called without additional arguments) can be
52##' used to inspect the linear constraints of a \code{lvm}-object.
53##'
54##' For backward compatibility the "$"-symbol can be used to fix parameters at
55##' a given value. E.g. to add a linear relationship between \code{y} and
56##' \code{x} with slope 2 to the model \code{m}, we can write
57##' \code{regression(m,"y") <- "x$2"}.  Similarily we can use the "@@"-symbol to
58##' name parameters. E.g. in a multiple regression we can force the parameters
59##' to be equal: \code{regression(m,"y") <- c("x1@@b","x2@@b")}.  Fixed parameters
60##' can be reset by fixing (with \$) them to \code{NA}.
61##'
62##' @aliases regression regression<- regression<-.lvm regression.lvm regfix
63##' regfix regfix<- regfix.lvm regfix<-.lvm
64##' @param object \code{lvm}-object.
65##' @param value A formula specifying the linear constraints or if
66##' \code{to=NULL} a \code{list} of parameter values.
67##' @param to Character vector of outcome(s) or formula object.
68##' @param from Character vector of predictor(s).
69##' @param fn Real function defining the functional form of predictors (for
70##' simulation only).
71##' @param messages Controls which messages are turned on/off (0: all off)
72##' @param additive If FALSE and predictor is categorical a non-additive effect is assumed
73##' @param y Alias for 'to'
74##' @param x Alias for 'from'
75##' @param quick Faster implementation without parameter constraints
76##' @param \dots Additional arguments to be passed to the low level functions
77##' @usage
78##' \method{regression}{lvm}(object = lvm(), to, from, fn = NA,
79##' messages = lava.options()$messages, additive=TRUE, y, x, value, ...)
80##' \method{regression}{lvm}(object, to=NULL, quick=FALSE, ...) <- value
81##' @return A \code{lvm}-object
82##' @note Variables will be added to the model if not already present.
83##' @author Klaus K. Holst
84##' @seealso \code{\link{intercept<-}}, \code{\link{covariance<-}},
85##' \code{\link{constrain<-}}, \code{\link{parameter<-}},
86##' \code{\link{latent<-}}, \code{\link{cancel<-}}, \code{\link{kill<-}}
87##' @keywords models regression
88##' @examples
89##'
90##' m <- lvm() ## Initialize empty lvm-object
91##' ### E(y1|z,v) = beta1*z + beta2*v
92##' regression(m) <- y1 ~ z + v
93##' ### E(y2|x,z,v) = beta*x + beta*z + 2*v + beta3*u
94##' regression(m) <- y2 ~ f(x,beta) + f(z,beta)  + f(v,2) + u
95##' ### Clear restriction on association between y and
96##' ### fix slope coefficient of u to beta
97##' regression(m, y2 ~ v+u) <- list(NA,"beta")
98##'
99##' regression(m) ## Examine current linear parameter constraints
100##'
101##' ## ## A multivariate model, E(yi|x1,x2) = beta[1i]*x1 + beta[2i]*x2:
102##' m2 <- lvm(c(y1,y2) ~ x1+x2)
103##'
104##' @export
105"regression<-" <- function(object,...,value) UseMethod("regression<-")
106
107##' @export
108regression.formula <- function(object,...) regression(lvm(),object,...)
109
110##' @export
111"regression<-.lvm" <- function(object, to=NULL, quick=FALSE, ..., value) {
112    dots <- list(...)
113    if (length(dots$additive)>0 && !dots$additive && !inherits(value,"formula")) {
114        regression(object,beta=value,...) <- to
115        return(object)
116    }
117    if (!is.null(to) || !is.null(dots$y)) {
118        regfix(object, to=to, ...) <- value
119        return(object)
120    } else  {
121        if (is.list(value)) {
122            for (v in value) {
123                regression(object,...) <- v
124            }
125            return(object)
126        }
127
128        if (inherits(value,"formula")) {
129            fff <- procformula(object,value,...)
130            object <- fff$object
131            lhs <- fff$lhs
132            xs <- fff$xs
133            ys <- fff$ys
134            res <- fff$res
135            X <- fff$X
136
137
138        if (fff$iscovar) {
139            ## return(covariance(object,var1=decomp.specials(lhs[[1]]),var2=X))
140            covariance(object) <- toformula(decomp.specials(lhs[[1]]),X)
141            return(object)
142        }
143        if (!is.null(lhs) && nchar(lhs[[1]])>2 && substr(lhs[[1]],1,2)=="v(") {
144            v <- update(value,paste(decomp.specials(lhs),"~."))
145            covariance(object,...) <- v
146            return(object)
147        }
148
149        if (length(lhs)==0) {
150            index(object) <- reindex(object)
151            return(object)
152        }
153
154        for (i in seq_len(length(ys))) {
155        y <- ys[i]
156        for (j in seq_len(length(xs))) {
157          if (length(res[[j]])>1) {
158            regfix(object, to=y[1], from=xs[j],...) <- res[[j]][2]
159          } else {
160            object <- regression(object,to=y[1],from=xs[j],...)
161          }
162        }
163      }
164      object$parpos <- NULL
165      return(object)
166    }
167    if (!is.list(value) | length(value)>2) stop("Value should contain names of outcome (to) and predictors (from)")
168    if (all(c("to","from")%in%names(value))) {
169
170      xval <- value$x; yval <- value$y
171    } else {
172      yval <- value[[1]]; xval <- value[[2]]
173    }
174    regression(object, to=yval, from=xval,...)
175  }
176}
177
178##' @export
179`regression` <-
180  function(object,to,from,...) UseMethod("regression")
181
182##' @export
183`regression.lvm` <-
184    function(object=lvm(),to,from,fn=NA,messages=lava.options()$messages,
185      additive=TRUE, y,x,value,...) {
186        if (!missing(y)) {
187            if (inherits(y,"formula")) y <- all.vars(y)
188            to <- y
189        }
190        if (!missing(x)) {
191            if (inherits(x,"formula")) x <- all.vars(x)
192            from <- x
193        }
194        if (!additive) {
195            if (!inherits(to,"formula")) to <- toformula(to,from)
196            x <- attributes(getoutcome(to))$x
197            K <- object$attributes$nordinal[x]
198            if (is.null(K) || is.na(K)) {
199                K <- list(...)$K
200                if (is.null(K)) stop("Supply number of categories, K (or use method 'categorical' before calling 'regression').")
201                object <- categorical(object,x,...)
202            }
203            dots <- list(...);
204            dots$K <- K
205            dots$x <- object
206            dots$formula <- to
207            dots$regr.only <- TRUE
208            object <- do.call("categorical",dots)
209            return(object)
210        }
211
212        if (missing(to)) {
213            return(regfix(object))
214        }
215        if (inherits(to,"formula")) {
216          if (!missing(value)) {
217                regression(object,to,messages=messages,...) <- value
218            } else {
219                regression(object,messages=messages,...) <- to
220            }
221            object$parpos <- NULL
222            return(object)
223        }
224        if (is.list(to)) {
225            for (t in to)
226                regression(object,messages=messages,...) <- t
227            object$parpos <- NULL
228            return(object)
229        }
230
231        sx <- strsplit(from,"@")
232        xx <- sapply(sx, FUN=function(i) i[1])
233        ps <- sapply(sx, FUN=function(i) i[2])
234        sx <- strsplit(xx,"$",fixed=TRUE)
235        xs <- sapply(sx, FUN=function(i) i[1])
236        fix <- char2num(sapply(sx, FUN=function(i) i[2]))
237        allv <- index(object)$vars
238
239        object <- addvar(object, c(to,xs), messages=messages, reindex=FALSE)
240
241        for (i in to)
242            for (j in xs) {
243                object$M[j,i] <- 1
244                if (!is.na(fn))
245                    functional(object,j,i) <- fn
246            }
247
248        if (lava.options()$exogenous) {
249            newexo <- setdiff(xs,c(to,allv))
250            exo <- exogenous(object)
251            if (length(newexo)>0)
252                exo <- unique(c(exo,newexo))
253            exogenous(object) <- setdiff(exo,to)
254        }
255
256        if (lava.options()$debug) {
257            print(object$fix)
258        }
259        object$fix[xs,to] <- fix
260        object$par[xs,to] <- ps
261        object$parpos <- NULL
262
263        index(object) <- reindex(object)
264        return(object)
265    }
266
267