##' Simulate model
##'
##' Simulate data from a general SEM model including non-linear effects and
##' general link and distribution of variables.
##'
##' @aliases sim sim.lvmfit sim.lvm
##' simulate.lvmfit simulate.lvm
##' transform<- transform<-.lvm transform.lvm
##' functional functional<- functional.lvm functional<-.lvm
##' distribution distribution distribution<- distribution.lvm distribution<-.lvm
##' heavytail heavytail<-
##' weibull.lvm
##' binomial.lvm
##' poisson.lvm
##' uniform.lvm
##' multinomial.lvm
##' beta.lvm
##' normal.lvm mvn.lvm
##' lognormal.lvm
##' gaussian.lvm
##' GM2.lvm
##' GM3.lvm
##' probit.lvm
##' logit.lvm
##' pareto.lvm
##' student.lvm
##' chisq.lvm
##' coxGompertz.lvm
##' coxWeibull.lvm
##' coxExponential.lvm
##' aalenExponential.lvm
##' Gamma.lvm gamma.lvm
##' loggamma.lvm
##' categorical categorical<-
##' threshold.lvm
##' ones.lvm Binary.lvm binary.lvm
##' Sequence.lvm
##' none.lvm constant.lvm id.lvm
##' @usage
##' \method{sim}{lvm}(x, n = NULL, p = NULL, normal = FALSE, cond = FALSE,
##' sigma = 1, rho = 0.5, X = NULL, unlink=FALSE, latent=TRUE,
##' use.labels = TRUE, seed=NULL, ...)
##' @param x Model object
##' @param n Number of simulated values/individuals
##' @param p Parameter value (optional)
##' @param normal Logical indicating whether to simulate data from a
##' multivariate normal distribution conditional on exogenous variables hence
##' ignoring functional/distribution definition
##' @param cond for internal use
##' @param sigma Default residual variance (1)
##' @param rho Default covariance parameter (0.5)
##' @param X Optional matrix of fixed values of variables (manipulation)
##' @param unlink Return Inverse link transformed data
##' @param latent Include latent variables (default TRUE)
##' @param use.labels convert categorical variables to factors before applying transformation
##' @param seed Random seed
##' @param \dots Additional arguments to be passed to the low level functions
##' @author Klaus K. Holst
##' @keywords models datagen regression
##' @export
##' @examples
##' ##################################################
##' ## Logistic regression
##' ##################################################
##' m <- lvm(y~x+z)
##' regression(m) <- x~z
##' distribution(m,~y+z) <- binomial.lvm("logit")
##' d <- sim(m,1e3)
##' head(d)
##'
##' e <- estimate(m,d,estimator="glm")
##' e
##' ## Simulate a few observation from estimated model
##' sim(e,n=5)
##'
##' ##################################################
##' ## Poisson
##' ##################################################
##' distribution(m,~y) <- poisson.lvm()
##' d <- sim(m,1e4,p=c(y=-1,"y~x"=2,z=1))
##' head(d)
##' estimate(m,d,estimator="glm")
##' mean(d$z); lava:::expit(1)
##'
##' summary(lm(y~x,sim(lvm(y[1:2]~4*x),1e3)))
##'
##' ##################################################
##' ### Gamma distribution
##' ##################################################
##' m <- lvm(y~x)
##' distribution(m,~y+x) <- list(Gamma.lvm(shape=2),binomial.lvm())
##' intercept(m,~y) <- 0.5
##' d <- sim(m,1e4)
##' summary(g <- glm(y~x,family=Gamma(),data=d))
##' \dontrun{MASS::gamma.shape(g)}
##'
##' args(lava::Gamma.lvm)
##' distribution(m,~y) <- Gamma.lvm(shape=2,log=TRUE)
##' sim(m,10,p=c(y=0.5))[,"y"]
##'
##' ##################################################
##' ### Beta
##' ##################################################
##' m <- lvm()
##' distribution(m,~y) <- beta.lvm(alpha=2,beta=1)
##' var(sim(m,100,"y,y"=2))
##' distribution(m,~y) <- beta.lvm(alpha=2,beta=1,scale=FALSE)
##' var(sim(m,100))
##'
##' ##################################################
##' ### Transform
##' ##################################################
##' m <- lvm()
##' transform(m,xz~x+z) <- function(x) x[1]*(x[2]>0)
##' regression(m) <- y~x+z+xz
##' d <- sim(m,1e3)
##' summary(lm(y~x+z + x*I(z>0),d))
##'
##' ##################################################
##' ### Non-random variables
##' ##################################################
##' m <- lvm()
##' distribution(m,~x+z+v+w) <- list(Sequence.lvm(0,5),## Seq. 0 to 5 by 1/n
##' Binary.lvm(), ## Vector of ones
##' Binary.lvm(0.5), ## 0.5n 0, 0.5n 1
##' Binary.lvm(interval=list(c(0.3,0.5),c(0.8,1))))
##' sim(m,10)
##'
##' ##################################################
##' ### Cox model
##' ### piecewise constant hazard
##' ################################################
##' m <- lvm(t~x)
##' rates <- c(1,0.5); cuts <- c(0,5)
##' ## Constant rate: 1 in [0,5), 0.5 in [5,Inf)
##' distribution(m,~t) <- coxExponential.lvm(rate=rates,timecut=cuts)
##'
##' \dontrun{
##' d <- sim(m,2e4,p=c("t~x"=0.1)); d$status <- TRUE
##' plot(timereg::aalen(survival::Surv(t,status)~x,data=d,
##' resample.iid=0,robust=0),spec=1)
##' L <- approxfun(c(cuts,max(d$t)),f=1,
##' cumsum(c(0,rates*diff(c(cuts,max(d$t))))),
##' method="linear")
##' curve(L,0,100,add=TRUE,col="blue")
##' }
##'
##' ##################################################
##' ### Cox model
##' ### piecewise constant hazard, gamma frailty
##' ##################################################
##' m <- lvm(y~x+z)
##' rates <- c(0.3,0.5); cuts <- c(0,5)
##' distribution(m,~y+z) <- list(coxExponential.lvm(rate=rates,timecut=cuts),
##' loggamma.lvm(rate=1,shape=1))
##' \dontrun{
##' d <- sim(m,2e4,p=c("y~x"=0,"y~z"=0)); d$status <- TRUE
##' plot(timereg::aalen(survival::Surv(y,status)~x,data=d,
##' resample.iid=0,robust=0),spec=1)
##' L <- approxfun(c(cuts,max(d$y)),f=1,
##' cumsum(c(0,rates*diff(c(cuts,max(d$y))))),
##' method="linear")
##' curve(L,0,100,add=TRUE,col="blue")
##' }
##' ## Equivalent via transform (here with Aalens additive hazard model)
##' m <- lvm(y~x)
##' distribution(m,~y) <- aalenExponential.lvm(rate=rates,timecut=cuts)
##' distribution(m,~z) <- Gamma.lvm(rate=1,shape=1)
##' transform(m,t~y+z) <- prod
##' sim(m,10)
##' ## Shared frailty
##' m <- lvm(c(t1,t2)~x+z)
##' rates <- c(1,0.5); cuts <- c(0,5)
##' distribution(m,~y) <- aalenExponential.lvm(rate=rates,timecut=cuts)
##' distribution(m,~z) <- loggamma.lvm(rate=1,shape=1)
##' \dontrun{
##' mets::fast.reshape(sim(m,100),varying="t")
##' }
##'
##' ##################################################
##' ### General multivariate distributions
##' ##################################################
##' \dontrun{
##' m <- lvm()
##' distribution(m,~y1+y2,oratio=4) <- VGAM::rbiplackcop
##' ksmooth2(sim(m,1e4),rgl=FALSE,theta=-20,phi=25)
##'
##' m <- lvm()
##' distribution(m,~z1+z2,"or1") <- VGAM::rbiplackcop
##' distribution(m,~y1+y2,"or2") <- VGAM::rbiplackcop
##' sim(m,10,p=c(or1=0.1,or2=4))
##' }
##'
##' m <- lvm()
##' distribution(m,~y1+y2+y3,TRUE) <- function(n,...) rmvn0(n,sigma=diag(3)+1)
##' var(sim(m,100))
##'
##' ## Syntax also useful for univariate generators, e.g.
##' m <- lvm(y~x+z)
##' distribution(m,~y,TRUE) <- function(n) rnorm(n,mean=1000)
##' sim(m,5)
##' distribution(m,~y,"m1",0) <- rnorm
##' sim(m,5)
##' sim(m,5,p=c(m1=100))
##'
##' ##################################################
##' ### Regression design in other parameters
##' ##################################################
##' ## Variance heterogeneity
##' m <- lvm(y~x)
##' distribution(m,~y) <- function(n,mean,x) rnorm(n,mean,exp(x)^.5)
##' if (interactive()) plot(y~x,sim(m,1e3))
##' ## Alternaively, calculate the standard error directly
##' addvar(m) <- ~sd ## If 'sd' should be part of the resulting data.frame
##' constrain(m,sd~x) <- function(x) exp(x)^.5
##' distribution(m,~y) <- function(n,mean,sd) rnorm(n,mean,sd)
##' if (interactive()) plot(y~x,sim(m,1e3))
##'
##' ## Regression on variance parameter
##' m <- lvm()
##' regression(m) <- y~x
##' regression(m) <- v~x
##' ##distribution(m,~v) <- 0 # No stochastic term
##' ## Alternative:
##' ## regression(m) <- v[NA:0]~x
##' distribution(m,~y) <- function(n,mean,v) rnorm(n,mean,exp(v)^.5)
##' if (interactive()) plot(y~x,sim(m,1e3))
##'
##' ## Regression on shape parameter in Weibull model
##' m <- lvm()
##' regression(m) <- y ~ z+v
##' regression(m) <- s ~ exp(0.6*x-0.5*z)
##' distribution(m,~x+z) <- binomial.lvm()
##' distribution(m,~cens) <- coxWeibull.lvm(scale=1)
##' distribution(m,~y) <- coxWeibull.lvm(scale=0.1,shape=~s)
##' eventTime(m) <- time ~ min(y=1,cens=0)
##'
##' if (interactive()) {
##' d <- sim(m,1e3)
##' require(survival)
##' (cc <- coxph(Surv(time,status)~v+strata(x,z),data=d))
##' plot(survfit(cc) ,col=1:4,mark.time=FALSE)
##' }
##'
##' ##################################################
##' ### Categorical predictor
##' ##################################################
##' m <- lvm()
##' ## categorical(m,K=3) <- "v"
##' categorical(m,labels=c("A","B","C")) <- "v"
##'
##' regression(m,additive=FALSE) <- y~v
##' \dontrun{
##' plot(y~v,sim(m,1000,p=c("y~v:2"=3)))
##' }
##'
##' m <- lvm()
##' categorical(m,labels=c("A","B","C"),p=c(0.5,0.3)) <- "v"
##' regression(m,additive=FALSE,beta=c(0,2,-1)) <- y~v
##' ## equivalent to:
##' ## regression(m,y~v,additive=FALSE) <- c(0,2,-1)
##' regression(m,additive=FALSE,beta=c(0,4,-1)) <- z~v
##' table(sim(m,1e4)$v)
##' glm(y~v, data=sim(m,1e4))
##' glm(y~v, data=sim(m,1e4,p=c("y~v:1"=3)))
##'
##' transform(m,v2~v) <- function(x) x=='A'
##' sim(m,10)
##'
##' ##################################################
##' ### Pre-calculate object
##' ##################################################
##' m <- lvm(y~x)
##' m2 <- sim(m,'y~x'=2)
##' sim(m,10,'y~x'=2)
##' sim(m2,10) ## Faster
"sim" <- function(x,...) UseMethod("sim")
##' @export
sim.lvmfit <- function(x,n=nrow(model.frame(x)),p=pars(x),xfix=TRUE,...) {
m <- Model(x)
if ((nrow(model.frame(x))==n) & xfix) {
X <- exogenous(x)
mydata <- model.frame(x)
for (pred in X) {
distribution(m, pred) <- list(mydata[,pred])
}
}
sim(m,n=n,p=p,...)
}
##' @export
sim.lvm <- function(x,n=NULL,p=NULL,normal=FALSE,cond=FALSE,sigma=1,rho=.5,
X=NULL,unlink=FALSE,latent=TRUE,use.labels=TRUE,seed=NULL,...) {
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
v.env <- c("A","M","P","PP","PPdiag","xx","vv","mdist","mdistnam","mii",
"nn","mu","xf","xfix","X",
"vartrans","multitrans","multitrans.idx",
"X.idx","ii.mvn","xconstrain.idx","xconstrain",
"xconstrain.par","covparnames","exo_constrainY")
setup <- is.null(n) && is.null(X) ## Save environment (variables v.env) and return sim object
loadconfig <- !is.null(x$sim.env) && !setup && (length(list(...))==0 && length(p)==0)
Yfix <- NULL
if (loadconfig) {
for (v in setdiff(v.env,"X")) assign(v, x$sim.env[[v]])
if (is.null(X)) X <- x$sim.env[['X']]
} else {
if (!is.null(n) && n<1) return(NULL)
p <- c(p,unlist(list(...)))
xx <- exogenous(x)
if (!is.null(X)) {
if (is.null(n))
n <- nrow(X)
if (!is.null(colnames(X))) {
yfix <- setdiff(colnames(X),xx)
if (length(yfix)>0) Yfix <- X[,yfix,drop=FALSE]
xx0 <- intersect(xx,colnames(X))
if (length(xx0)>0)
X <- as.matrix(X[,xx0,drop=FALSE])
}
} else {
if (!is.null(p)) {
i1 <- unique(na.omit(c(match(names(p),xx),
match(names(p),paste0(xx,lava.options()$symbol[2],xx)))))
covariance(x) <- xx[i1]
}
}
## index(x) <- reindex(x)
vv <- vars(x)
nn <- setdiff(vv,parameter(x))
mu <- unlist(lapply(x$mean, function(l) ifelse(is.na(l)|is.character(l),0,l)))
xf <- intersect(unique(parlabels(x)),xx)
xfix <- c(randomslope(x),xf); if (length(xfix)>0) normal <- FALSE
## Match parameter names
if ((!is.null(names(p)) && all(!is.na(names(p)))) || length(p)!=(index(x)$npar+index(x)$npar.mean+index(x)$npar.ex) | is.null(names(p))) {
nullp <- is.null(p)
p0 <- p
ep <- NULL
ei <- which(index(x)$e1==1)
if (length(ei)>0)
ep <- unlist(x$expar)[ei]
p <- c(rep(1, index(x)$npar+index(x)$npar.mean),ep)
p[seq_len(index(x)$npar.mean)] <- 0
p[index(x)$npar.mean + variances(x)] <- sigma
p[index(x)$npar.mean + offdiags(x)] <- rho
if (!nullp) {
c1 <- coef(x,mean=TRUE,fix=FALSE)
c2 <- coef(x,mean=TRUE,fix=FALSE,labels=TRUE)
idx1 <- na.omit(match(names(p0),c1))
idx11 <- na.omit(match(names(p0),c2))
idx2 <- na.omit(which(names(p0)%in%c1))
idx22 <- na.omit(which(names(p0)%in%c2))
if (length(idx1)>0 && !any(is.na(idx1)))
p[idx1] <- p0[idx2]
if (length(idx11)>0 && !any(is.na(idx11)))
p[idx11] <- p0[idx22]
}
}
M <- modelVar(x,p,data=NULL)
A <- M$A; P <- M$P
if (!is.null(M$v)) mu <- M$v
## Square root of residual variance matrix
PP <- with(svd(P), v%*%diag(sqrt(d),nrow=length(d))%*%t(u))
## Multivariate distributions
mdist <- distribution(x,multivariate=TRUE)$var
mdistnam <- names(mdist)
mii <- match(mdistnam,vars(x))
if (length(distribution(x))>0 ) {
ii <- match(names(distribution(x)),vv)
ii.mvn <- setdiff(seq(ncol(P)),c(ii,mii))
} else {
ii.mvn <- seq(ncol(P))
}
PPdiag <- sum(abs(offdiag(PP[ii.mvn,ii.mvn,drop=FALSE])^2))<1e-20
}
if (!setup) {
E <- matrix(0,ncol=ncol(P),nrow=n)
if (length(ii.mvn)>0) {
## Error term for conditional normal distributed variables
if (PPdiag) {
for (i in ii.mvn) E[,i] <- rnorm(n,sd=PP[i,i])
} else {
E[,ii.mvn] <- matrix(rnorm(length(ii.mvn)*n),ncol=length(ii.mvn))%*%PP[ii.mvn,ii.mvn,drop=FALSE]
}
}
if (length(mdistnam)>0) {
fun <- distribution(x,multivariate=TRUE)$fun
for (i in seq_along(fun)) {
mv <- names(which(unlist(mdist)==i))
ii <- match(mv,vv)
E[,ii] <- distribution(x,multivariate=TRUE)$fun[[i]](n,p=p,object=x) # ,...)
}
}
## Simulate exogenous variables (covariates)
res <- matrix(0,ncol=length(nn),nrow=n)
colnames(res) <- nn
}
if (!loadconfig) {
vartrans <- names(x$attributes$transform)
multitrans <- multitrans.idx <- NULL
if (length(x$attributes$multitransform)>0) {
multitrans <- unlist(lapply(x$attributes$multitransform,function(z) z$y))
for (i in (seq_along(x$attributes$multitransform))) {
multitrans.idx <- c(multitrans.idx,rep(i,length(x$attributes$multitransform[[i]]$y)))
}
}
xx <- unique(c(exogenous(x, latent=FALSE, index=TRUE),xfix))
xx <- setdiff(xx,vartrans)
X.idx <- match(xx,vv)
}
if (!setup) {
res[,X.idx] <- t(mu[X.idx]+t(E[,X.idx]))
if (is.null(X) || NCOL(X)<length(xx)) {
if (!is.null(xx) && length(xx)>0)
for (i in seq_along(xx)) {
mu.x <- mu[X.idx[i]]
dist.x <- distribution(x,xx[i])[[1]]
if (is.list(dist.x) && is.function(dist.x[[1]])) dist.x <- dist.x[[1]]
if (is.list(dist.x)) {
dist.x <- dist.x[[1]]
if (length(dist.x)==1) dist.x <- rep(dist.x,n)
}
if (is.function(dist.x)) {
res[,X.idx[i]] <- dist.x(n=n,mu=mu.x,var=P[X.idx[i],X.idx[i]])
} else {
if (is.null(dist.x) || any(is.na(dist.x))) {
} else {
if (length(dist.x)!=n) stop("'",vv[X.idx[i]], "' fixed at length ", length(dist.x)," != ",n,".")
res[,X.idx[i]] <- dist.x ## Deterministic
}
}
}
}
if (!is.null(X)) {
ii <- match(colnames(X),vv)
res0 <- res
for (i in seq(ncol(X)))
res[,ii[i]] <- X[,i]
}
}
simuled <- c(xx)
resunlink <- NULL
if (unlink) {
resunlink <- res
}
if ( normal | ( is.null(distribution(x)) & is.null(functional(x)) & is.null(constrain(x))) ) {
if(cond) { ## Simulate from conditional distribution of Y given X
mypar <- pars(x,A,P,mu)
Ey.x <- predict(x, mypar, data.frame(res))
Vy.x <- attributes(Ey.x)$cond.var
PP <- with(svd(Vy.x), v%*%diag(sqrt(d),nrow=length(d))%*%t(u))
yy <- Ey.x + matrix(n*ncol(Vy.x),ncol=ncol(Vy.x))%*%PP
res <- cbind(yy, res[,xx]); colnames(res) <- c(colnames(Vy.x),xx)
return(res)
}
## Simulate from sim. distribution (Y,X) (mv-normal)
I <- diag(nrow=length(nn))
IAi <- Inverse(I-t(A))
colnames(E) <- vv
dd <- t(apply(heavytail_sim_hook(x,E),1,function(x) x+mu))
res <- dd%*%t(IAi)
colnames(res) <- vv
} else {
if (!loadconfig) {
xc <- index(x)$vars
xconstrain.idx <- unlist(lapply(lapply(constrain(x),function(z) attributes(z)$args),function(z) length(intersect(z,xc))>0))
xconstrain <- intersect(unlist(lapply(constrain(x),function(z) attributes(z)$args)),xc)
xconstrain.par <- names(xconstrain.idx)[xconstrain.idx]
covparnames <- unique(as.vector(covariance(x)$labels))
exo_constrainY <- intersect(exogenous(x),names(x$constrainY))
}
if (setup) {
sim.env <- c()
sim.env[v.env] <- list(NULL)
for (v in v.env) if (!is.null(get(v))) sim.env[[v]] <- get(v)
x$sim.env <- sim.env
return(x)
}
if (length(xconstrain)>0) {
for (i in which(xconstrain.idx)) {
if (!(names(xconstrain.idx[i]) %in% nn)) { ## A parameter and not a variable
ff <- constrain(x)[[i]]
myargs <- attributes(ff)$args
D <- matrix(0,n,length(myargs))
for (j in seq_len(ncol(D))) {
if (myargs[j]%in%xconstrain)
D[,j] <- res[,myargs[j]]
else
D[,j] <- M$parval[[myargs[j]]]
}
val <- try(val <- ff(D), silent=TRUE)
if (inherits(val,"try-error") || NROW(val)<n) apply(D,1,ff)
res <- cbind(res, val); colnames(res)[ncol(res)] <- names(xconstrain.idx)[i]
}
}
}
if (any(xconstrain.par%in%covparnames)) {
mu0 <- rep(0,ncol(P))
P0 <- P
E <- t(sapply(seq_len(n),function(idx) {
for (i in intersect(xconstrain.par,covparnames)) {
P0[covariance(x)$labels==i] <- res[idx,i]
}
PP <- with(svd(P0), v%*%diag(sqrt(d),nrow=length(d))%*%t(u))
return(mu0+rbind(rnorm(ncol(P0)))%*%PP)
}))
}
colnames(E) <- vv
E <- heavytail_sim_hook(x,E)
## Non-linear regression components
xconstrain <- c()
for (i in seq_len(length(constrain(x)))) {
z <- constrain(x)[[i]]
xx <- intersect(attributes(z)$args,manifest(x))
if (length(xx)>0) {
warg <- setdiff(attributes(z)$args,xx)
wargidx <- which(attributes(z)$args%in%warg)
exoidx <- which(attributes(z)$args%in%xx)
parname <- names(constrain(x))[i]
y <- names(which(unlist(lapply(intercept(x),function(x) x==parname))))
el <- list(i,y,parname,xx,exoidx,warg,wargidx,z)
names(el) <- c("idx","endo","parname","exo","exoidx","warg","wargidx","func")
xconstrain <- c(xconstrain,list(el))
}
}
yconstrain <- unlist(lapply(xconstrain,function(x) x$endo))
for (i in exo_constrainY) {
cc <- x$constrainY[[i]]
args <- cc$args
args <- if (is.null(args) || length(args)==0) res[,i] else res[,args]
res[,i] <- cc$fun(args,p) # ,...)
}
res <- data.frame(res)
if (length(vartrans)>0) {
parvals <- parpos(x)$parval
parvalsnam <- setdiff(names(parvals),xx)
if (length(parvalsnam)>0) {
Parvals <- p[unlist(parvals)];
res <- cbind(res,
cbind(rep(1,nrow(res)))%x%rbind(Parvals))
colnames(res)[seq(length(Parvals))+ncol(res)-length(Parvals)] <-
names(parvals)
}
}
leftovers <- c()
itercount <- 0
while (length(simuled)<length(nn)) {
leftoversPrev <- leftovers
leftovers <- setdiff(nn,simuled)
if (!is.null(leftoversPrev) && length(leftoversPrev)==length(leftovers)) {
if (itercount>0)
stop("Infinite loop (feedback).")
itercount <- itercount+1
}
for (i in leftovers) {
if (length(Yfix)>0 && i %in% colnames(Yfix)) {
if (NROW(Yfix) == 1) {
res[,i] <- rep(Yfix[,i],length.out=nrow(res))
} else
res[,i] <- rep(Yfix[,i,drop=TRUE],length.out=nrow(res))
simuled <- c(simuled,i)
next
}
if (i%in%vartrans) {
xtrans <- x$attributes$transform[[i]]$x
if (all(xtrans%in%c(simuled,names(parvals)))) {
xtr <- res[,xtrans,drop=FALSE]
if (use.labels) {
lb <- x$attributes$labels
lb.idx <- na.omit(match(names(lb),xtrans))
## For categorical variables turn them into factors so we can
## use the actual labels in function calls/transform
if (length(lb.idx)>0) {
xtr <- as.data.frame(xtr)
for (lb0 in lb.idx) {
lab <- lb[[names(xtr)[lb0]]]
xtr[,lb0] <- factor(xtr[,lb0],levels=seq_along(lab)-1,labels=lab)
}
}
}
suppressWarnings(yy <- with(x$attributes$transform[[i]], fun(xtr))) ##fun(res[,xtrans])))
if (NROW(yy) != NROW(res)) { ## apply row-wise
res[,i] <- with(x$attributes$transform[[i]], ##apply(res[,xtrans,drop=FALSE],1,fun))
apply(xtr,1,fun))
} else {
colnames(yy) <- NULL
res[,i] <- yy
}
simuled <- c(simuled,i)
}
} else if (i%in%multitrans) {
idx0 <- match(i,multitrans)
idx <- multitrans.idx[idx0]
mtval <- x$attributes$multitransform[[idx]]
if (all(mtval$x%in%simuled)) {
res[,mtval$y] <- mtval$fun(res[,mtval$x])
simuled <- c(simuled,mtval$y)
break;
}
} else {
ipos <- which(yconstrain%in%i)
if (length(ipos)==0 || all(xconstrain[[ipos]]$exo%in%simuled)) {
pos <- match(i,vv)
relations <- colnames(A)[A[,pos]!=0]
simvars <- x$attributes$simvar[[i]]
dist.i <- distribution(x,i)[[1]] ## User-specified distribution function
dist.xx <- NULL
if (is.function(dist.i)) {
dist.args0 <- names(formals(dist.i))
dist.args <- setdiff(dist.args0,c("n","mean","mu","var","..."))
dist.xx <- intersect(names(res),dist.args) ## Variables influencing distribution
}
if (all(c(relations,simvars,dist.xx)%in%simuled)) { ## Only depending on already simulated variables
if (x$mean[[pos]]%in%xconstrain.par && length(ipos)==0) {
mu.i <- res[,x$mean[[pos]] ]
} else {
mu.i <- mu[pos]
}
if (length(ipos)>0) {
pp <- unlist(M$parval[xconstrain[[ipos]]$warg])
myidx <- with(xconstrain[[ipos]],order(c(wargidx,exoidx)))
## myidx <- with(xconstrain[[ipos]],
## match(attr(func,"args"), c(warg,exo)))
X <- with(xconstrain[[ipos]],
if (length(pp)>0)
cbind(rbind(pp)%x%cbind(rep(1,nrow(res))),
res[,exo,drop=FALSE])
else res[,exo,drop=FALSE])
yy <- try(with(xconstrain[[ipos]],
func(X[,myidx,drop=FALSE])),silent=TRUE)
if (NROW(yy) != NROW(res)) { ## apply row-wise
mu.i <- #mu.i +
with(xconstrain[[ipos]],
apply(res[,exo,drop=FALSE],1,
function(x) func(
unlist(c(pp,x))[myidx])))
} else {
mu.i <- ##mu.i+
yy
}
}
for (From in relations) {
f <- functional(x,i,From)[[1]]
if (!is.function(f))
f <- function(x,...) x
reglab <- regfix(x)$labels[From,pos]
if (reglab%in%c(xfix,xconstrain.par)) {
if (is.function(f)) {
if (length(formals(f))>1) {
mu.i <- mu.i + res[,reglab]*f(res[,From],p)
} else {
mu.i <- mu.i + res[,reglab]*f(res[,From])
}
} else mu.i <- mu.i + res[,reglab]*res[,From]
}
else {
if (is.function(f)) {
if (length(formals(f))>1) {
mu.i <- mu.i + A[From,pos]*f(res[,From],p)
} else {
mu.i <- mu.i + A[From,pos]*f(res[,From])
}
} else mu.i <- mu.i + A[From,pos]*res[,From]
}
}
if (!is.function(dist.i)) {
res[,pos] <- mu.i + E[,pos]
if (unlink)
resunlink[,pos] <- res[,pos]
}
else {
if (length(simvars)>0) { ## Depends on mu and also on other variables (e.g. time-depending effect)
if (length(mu.i)==1) mu.i <- rep(mu.i,n)
mu.i <- cbind("m0"=mu.i,res[,simvars,drop=FALSE])
}
new.args <- list(n=n)
mu.arg <- intersect(c("mean","mu"),dist.args0)
if (length(mu.arg)>0) {
new.args <- c(new.args,list(mu.i))
names(new.args)[length(new.args)] <- mu.arg[1]
}
var.arg <- intersect(c("var"),dist.args0)
if (length(var.arg)>0) {
new.args <- c(new.args,list(P[pos,pos]))
names(new.args)[length(new.args)] <- var.arg[1]
}
for (jj in dist.xx) {
new.args <- c(new.args,list(res[,jj,drop=TRUE]))
names(new.args)[length(new.args)] <- jj
}
res[,pos] <- do.call(dist.i,new.args)
if (unlink)
resunlink[,pos] <- mu.i
}
if (length(x$constrainY)>0 && i%in%names(x$constrainY)) {
cc <- x$constrainY[[i]]
args <- cc$args
args <- if (length(args)==0)
res[,pos]
else {
ii <- intersect(names(M$parval),args)
args0 <- args
args <- res[,intersect(args0,colnames(res)),drop=FALSE]
if (length(ii)>0) {
pp <- rbind(unlist(M$parval[ii]))%x%cbind(rep(1,n))
colnames(pp) <- ii
args <- cbind(res,pp)[,args0,drop=FALSE]
}
args
}
res[,pos] <- cc$fun(args,p) # ,...)
}
simuled <- c(simuled,i)
}
}
}
}
}
res <- res[,nn,drop=FALSE]
}
res <- as.data.frame(res)
myhooks <- gethook("sim.hooks")
for (f in myhooks) {
res <- do.call(f, list(x=x,data=res,p=p,modelpar=M))
}
if (unlink) res <- resunlink
res <- as.data.frame(res)
self <- x$attributes$selftransform
for (v in names(self)) {
res[,v] <- self[[v]](res[,v])
}
if (!latent && length(latent(x))>0) return(subset(res[,-which(colnames(res)%in%latent(x))]))
return(res)
}
##' @export
simulate.lvm <- function(object,nsim,seed=NULL,...) {
sim(object,nsim,seed=seed,...)
}
##' @export
simulate.lvmfit <- function(object,nsim,seed=NULL,...) {
sim(object,nsim,seed=seed,...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.