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