1##' Wrapper function for mclapply
2##'
3##' @export
4##' @param x function or 'sim' object
5##' @param R Number of replications or data.frame with parameters
6##' @param f Optional function (i.e., if x is a matrix)
7##' @param colnames Optional column names
8##' @param messages Messages
9##' @param seed (optional) Seed (needed with cl=TRUE)
10##' @param args (optional) list of named arguments passed to (mc)mapply
11##' @param iter If TRUE the iteration number is passed as first argument to (mc)mapply
12##' @param ... Additional arguments to (mc)mapply
13##' @aliases sim.default as.sim
14##' @seealso summary.sim plot.sim print.sim
15##' @examples
16##' m <- lvm(y~x+e)
17##' distribution(m,~y) <- 0
18##' distribution(m,~x) <- uniform.lvm(a=-1.1,b=1.1)
19##' transform(m,e~x) <- function(x) (1*x^4)*rnorm(length(x),sd=1)
20##'
21##' onerun <- function(iter=NULL,...,n=2e3,b0=1,idx=2) {
22##'     d <- sim(m,n,p=c("y~x"=b0))
23##'     l <- lm(y~x,d)
24##'     res <- c(coef(summary(l))[idx,1:2],
25##'              confint(l)[idx,],
26##'              estimate(l,only.coef=TRUE)[idx,2:4])
27##'     names(res) <- c("Estimate","Model.se","Model.lo","Model.hi",
28##'                     "Sandwich.se","Sandwich.lo","Sandwich.hi")
29##'     res
30##' }
31##' val <- sim(onerun,R=10,b0=1,messages=0)
32##' val
33##'
34##' val <- sim(val,R=40,b0=1) ## append results
35##' summary(val,estimate=c(1,1),confint=c(3,4,6,7),true=c(1,1))
36##'
37##' summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
38##' summary(val,estimate=c(1,1),se=c(2,5),true=c(1,1),names=c("Model","Sandwich"),confint=TRUE)
39##'
40##' if (interactive()) {
41##'     plot(val,estimate=1,c(2,5),true=1,names=c("Model","Sandwich"),polygon=FALSE)
42##'     plot(val,estimate=c(1,1),se=c(2,5),main=NULL,
43##'          true=c(1,1),names=c("Model","Sandwich"),
44##'          line.lwd=1,col=c("gray20","gray60"),
45##'          rug=FALSE)
46##'     plot(val,estimate=c(1,1),se=c(2,5),true=c(1,1),
47##'          names=c("Model","Sandwich"))
48##' }
49##''
50##' f <- function(a=1, b=1) {
51##'   rep(a*b, 5)
52##' }
53##' R <- Expand(a=1:3, b=1:3)
54##' sim(f, R)
55##' sim(function(a,b) f(a,b), 3, args=c(a=5,b=5))
56##' sim(function(iter=1,a=5,b=5) iter*f(a,b), iter=TRUE, R=5)
57sim.default <- function(x=NULL, R=100, f=NULL, colnames=NULL,
58                messages=lava.options()$messages,
59                seed=NULL, args=list(),
60                iter=FALSE, ...) {
61    stm <- proc.time()
62    oldtm <- rep(0,5)
63    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
64        runif(1)
65    if (is.null(seed))
66        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
67    else {
68        R.seed <- get(".Random.seed", envir = .GlobalEnv)
69        set.seed(seed)
70        RNGstate <- structure(seed, kind = as.list(RNGkind()))
71        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
72    }
73    olddata <- NULL
74    dots <- list(...)
75    mycall <- match.call(expand.dots=FALSE)
76    if (inherits(x, c("data.frame", "matrix"))) olddata <- x
77    if (inherits(x, "sim")) {
78        oldtm <- attr(x, "time")
79        oldcall <- attr(x, "call")
80        x <- attr(x, "f")
81        if (!is.null(f)) x <- f
82        ex <- oldcall[["..."]]
83        for (nn in setdiff(names(ex),names(dots))) {
84            dots[[nn]] <- ex[[nn]]
85            val <- list(ex[[nn]]); names(val) <- nn
86            mycall[["..."]] <- c(mycall[["..."]], list(val))
87        }
88
89    } else {
90        if (!is.null(f)) x <- f
91        if (!is.function(x)) stop("Expected a function or 'sim' object.")
92    }
93    if (is.null(x)) stop("Must give new function argument 'f'.")
94    res <- val <- NULL
95    on.exit({
96        if (is.null(colnames) && !is.null(val)) {
97            if (is.matrix(val[[1]])) {
98                colnames <- base::colnames(val[[1]])
99            } else {
100                colnames <- names(val[[1]])
101            }
102        }
103        base::colnames(res) <- colnames
104        if (!is.null(olddata)) res <- rbind(olddata,res)
105        attr(res, "call") <- mycall
106        attr(res, "f") <- x
107        class(res) <- c("sim", "matrix")
108        attr(res, "time") <- proc.time()-stm+oldtm
109        return(res)
110    })
111    parval_provided <- FALSE
112    if (inherits(R, c("matrix", "data.frame")) || length(R)>1) {
113        parval_provided <- TRUE
114        parval <- as.data.frame(R)
115        if (is.vector(R)) names(parval) <- NULL
116        else if (inherits(R,c("matrix", "data.frame")))
117          names(parval) <- colnames(R)
118        R <- NROW(parval)
119    } else {
120        parval <- as.data.frame(1:R)
121        names(parval) <- NULL
122    }
123
124    pb <- progressr::progressor(steps = R)
125    robx <- function(iter__, ...) {
126      pb()
127      tryCatch(x(...), error=function(e) NA)
128    }
129    if (iter) formals(robx)[[1]] <- NULL
130
131    pp <- c(as.list(parval),dots,
132            list(FUN=robx, SIMPLIFY=FALSE, MoreArgs=as.list(args), future.seed=TRUE))
133    val <- do.call(future.apply::future_mapply, pp)
134    res <- do.call(rbind, val)
135    if (is.null(res)) {
136      res <- matrix(NA, ncol=length(val[[1]]), nrow=R)
137    }
138    res
139}
140
141##' @export
142"[.sim" <- function (x, i, j, drop = FALSE) {
143    atr <- attributes(x)
144    if (!is.null(dim(x))) {
145        class(x) <- "matrix"
146    } else {
147        class(x) <- class(x)[-1]
148    }
149    x <- NextMethod("[",drop=drop)
150    atr.keep <- c("call","time")
151    if (missing(j)) atr.keep <- c(atr.keep,"f")
152    attributes(x)[atr.keep] <- atr[atr.keep]
153    if (!drop) class(x) <- c("sim",class(x))
154    x
155}
156
157##' @export
158"as.sim" <- function (object, name, ...) {
159    if (is.vector(object)) {
160        object <- (structure(cbind(object), class=c("sim", "matrix")))
161        if (!missing(name)) colnames(object) <- name
162        return(object)
163    }
164    structure(object, class=c("sim", class(object)))
165}
166
167Time <- function(sec,print=FALSE,...) {
168    h <- sec%/%3600
169    m0 <- (sec%%3600)
170    m <- m0%/%60
171    s <- m0%%60
172    res <- c(h=h,m=m,s=s)
173    if (print) {
174        if (h>0) cat(h,"h ",sep="")
175        if (m>0) cat(m,"m ",sep="")
176        cat(s,"s",sep="")
177        return(invisible(res))
178    }
179    return(res)
180}
181
182Print <- function(x,n=5,digits=max(3,getOption("digits")-3),...) {
183    mat <- !is.null(dim(x))
184    if (!mat) {
185        x <- cbind(x)
186        colnames(x) <- ""
187    }
188    if (is.null(rownames(x))) {
189        rownames(x) <- seq(nrow(x))
190    }
191    sep <- rbind("---"=rep('',ncol(x)))
192    if (n<1) {
193        print(x,quote=FALSE,digits=digits,...)
194    } else {
195        ## hd <- base::as.matrix(base::format(utils::head(x,n),digits=digits,...))
196        ## tl <- base::as.matrix(base::format(utils::tail(x,n),digits=digits,...))
197        ## print(rbind(hd,sep,tl),quote=FALSE,...)
198        if (NROW(x)<=(2*n)) {
199            hd <- base::format(utils::head(x,2*n),digits=digits,...)
200            print(hd, quote=FALSE,...)
201        } else {
202            hd <- base::format(utils::head(x,n),digits=digits,...)
203            tl <- base::format(utils::tail(x,n),digits=digits,...)
204            print(rbind(base::as.matrix(hd),sep,base::as.matrix(tl)),
205                  quote=FALSE,...)
206        }
207    }
208    invisible(x)
209}
210
211##' @export
212print.sim <- function(x,...) {
213    s <- summary(x,minimal=TRUE,...)
214    attr(x, "f") <- attr(x, "call") <- NULL
215    if (!is.null(dim(x))) {
216        class(x) <- "matrix"
217    }
218    Print(x,...)
219    cat("\n")
220    print(s,extra=FALSE,...)
221    return(invisible(x))
222}
223
224
225
226##' @export
227print.summary.sim <- function(x,group=list(c("^mean$","^sd$","^se$","^se/sd$","^coverage"),
228                                   c("^min$","^[0-9.]+%$","^max$"),
229                                   c("^na$","^missing$"),
230                                   c("^true$","^bias$","^rmse$")),
231                      lower.case=TRUE,
232                      na.print="",
233                      digits = max(3, getOption("digits") - 2),
234                      quote=FALSE,
235                      time=TRUE,
236                      extra=TRUE,
237                      ...) {
238    if (extra) {
239        cat(attr(x,"n")," replications",sep="")
240        if (time && !is.null(attr(x,"time"))) {
241            cat("\t\t\t\t\tTime: ")
242            Time(attr(x,"time")["elapsed"],print=TRUE)
243        }
244        cat("\n\n")
245    }
246
247    nn <- rownames(x)
248    if (lower.case)  nn <- tolower(nn)
249    gg <- lapply(group,
250                 function(x) unlist(lapply(x,function(v) grep(v,nn))))
251    gg <- c(gg,list(setdiff(seq_along(nn),unlist(gg))))
252
253    x0 <- c()
254    ng <- length(gg)
255    for (i in seq(ng)) {
256        x0 <- rbind(x0, x[gg[[i]],,drop=FALSE],
257        { if(i<ng && length(gg[[i+1]])>0) NA})
258    }
259
260    print(structure(x0,class="matrix")[,,drop=FALSE],digits=digits,quote=quote,na.print=na.print,...)
261    if (extra) cat("\n")
262    invisible(x)
263}
264
265
266##' Summary method for 'sim' objects
267##'
268##' Summary method for 'sim' objects
269##' @export
270##' @export summary.sim
271##' @param object sim object
272##' @param estimate (optional) columns with estimates
273##' @param se (optional) columns with standard error estimates
274##' @param confint (optional) list of pairs of columns with confidence limits
275##' @param true (optional) vector of true parameter values
276##' @param fun (optional) summary function
277##' @param names (optional) names of
278##' @param unique.names if TRUE, unique.names will be applied to column names
279##' @param minimal if TRUE, minimal summary will be returned
280##' @param level confidence level
281##' @param quantiles quantiles
282##' @param ... additional levels to lower-level functions
283summary.sim <- function(object,estimate=NULL,se=NULL,
284                confint=!is.null(se)&&!is.null(true),true=NULL,
285                fun,names=NULL,unique.names=TRUE,minimal=FALSE,
286                level=0.95,quantiles=c(0,.025,0.5,.975,1),...) {
287    if (is.list(estimate)) {
288        est <- estimate
289        if (is.null(names)) names <- base::names(est)
290        estimate <- c()
291        nse  <- is.null(se)
292        ntrue <- is.null(true)
293        elen <- unlist(lapply(est,length))
294        est <- lapply(est, function(e) c(e, rep(NA,max(elen)-length(e))))
295        for (e in est) {
296            estimate <- c(estimate,e[1])
297            if (length(e)>1 && nse) se <- c(se,e[2])
298            if (length(e)>2 && ntrue) true <- c(true,e[3])
299        }
300        cl <- match.call()
301        cl[c("estimate","se","true","names")] <- list(estimate,se,true,names)
302    }
303    if (minimal) {
304        fun <- function(x,se,confint,...) {
305            res <- c(Mean=mean(x,na.rm=TRUE),
306                    SD=sd(x,na.rm=TRUE))
307            if (!missing(se) && !is.null(se)) {
308                res <- c(res, c(SE=mean(se,na.rm=TRUE)))
309                res <- c(res, c("SE/SD"=res[["SE"]]/res[["SD"]]))
310            }
311            return(res)
312        }
313    }
314    mfun <- function(x,...) {
315        res <- c(mean(x,na.rm=TRUE),
316                 sd(x,na.rm=TRUE),
317                 if (length(quantiles)>0) quantile(x,quantiles,na.rm=TRUE),
318                 mean(is.na(x)))
319        if (length(quantiles)>0) {
320            nq <- paste0(quantiles*100,"%")
321            idx <- which(quantiles==1)
322            if (length(idx)>0) nq[idx] <- "Max"
323            idx <- which(quantiles==0)
324            if (length(idx)>0) nq[idx] <- "Min"
325        }
326        names(res) <- c("Mean","SD",
327                        if (length(quantiles)>0) nq,
328                        "Missing")
329        res
330    }
331    tm <- attr(object,"time")
332    N <- max(length(estimate),length(se),length(true))
333    if (!is.null(estimate)) estimate <- rep(estimate,length.out=N)
334    if (!is.null(se)) se <- rep(se,length.out=N)
335    if (!is.null(true)) {
336        if (is.null(estimate)) N <- ncol(object)
337        true <- rep(true,length.out=N)
338    }
339
340    if (!is.null(estimate) && is.character(estimate)) {
341        estimate <- match(estimate,colnames(object))
342    }
343    if (!missing(fun)) {
344        if (!is.null(estimate)) m.est <- object[,estimate,drop=FALSE]
345        else m.est <- object
346        m.se <- NULL
347        if (!is.null(se)) m.se <- object[,se,drop=FALSE]
348        m.ci <- NULL
349        if (!is.null(confint)) m.ci <- object[,confint,drop=FALSE]
350        res <- lapply(seq(ncol(m.est)),
351                      function(i,...) fun(m.est[,i,drop=TRUE],se=m.se[,i,drop=TRUE],confint=m.ci[,1:2+(i-1)*2],...,INDEX=i),...)
352        res <- matrix(unlist(res),nrow=length(res[[1]]),byrow=FALSE)
353        if (is.null(dim(res))) {
354            res <- rbind(res)
355        }
356        if (is.null(rownames(res))) {
357            rownames(res) <- names(fun(0,m.se,m.ci,INDEX=1,...))
358            if (is.null(rownames(res))) rownames(res) <- rep("",nrow(res))
359        }
360        if (is.null(colnames(res))) {
361            colnames(res) <- colnames(m.est)
362        }
363        return(structure(res,
364                    n=NROW(object),
365                    time=tm,
366                    class=c("summary.sim","matrix")))
367    }
368
369    if (!is.null(estimate)) {
370        est <- apply(object[,estimate,drop=FALSE],2,mfun)
371    } else {
372        est <- apply(object,2,mfun)
373    }
374
375    if (!is.null(true)) {
376        if (length(true)!=ncol(est)) {
377            ##stop("'true' should be of same length as 'estimate'.")
378            true <- rep(true,length.out=ncol(estimate))
379        }
380        est <- rbind(est,
381                     rbind(True=true),rbind(Bias=est["Mean",]-true),
382                     rbind(RMSE=((est["Mean",]-true)^2+(est["SD",])^2)^.5)
383                     )
384    }
385    if (!is.null(se)) {
386        if (is.character(se)) {
387            se <- match(se,colnames(object))
388        }
389        if (length(se)!=ncol(est)) stop("'se' should be of same length as 'estimate'.")
390        est <- rbind(est, SE=apply(object[,se,drop=FALSE],2,
391                                   function(x) val <- c(mean(x,na.rm=TRUE))))
392        est <- rbind(est,"SE/SD"=est["SE",]/est["SD",])
393    }
394    if (!is.null(confint) && (length(confint)>1 || confint)) {
395        if (is.character(confint)) {
396            confint <- match(confint,colnames(object))
397        }
398        if (length(confint)==1 && confint) {
399            if (is.null(se)) stop("Supply confidence limits or SE")
400            confint <- c()
401            pos <- ncol(object)
402            for (i in seq_along(estimate)) {
403                z <- 1-(1-level)/2
404                CI <- cbind(object[,estimate[i]]-qnorm(z)*object[,se[i]],
405                            object[,estimate[i]]+qnorm(z)*object[,se[i]])
406                colnames(CI) <- NULL
407                object <- cbind(object,CI)
408                confint <- c(confint,pos+1:2)
409                pos <- pos+2
410            }
411        }
412        if (length(confint)!=2*length(estimate)) stop("'confint' should be of length 2*length(estimate).")
413        Coverage <- c()
414        for (i in seq_along(estimate)) {
415            Coverage <- c(Coverage,
416                          mean((object[,confint[2*(i-1)+1]]<true[i]) & (object[,confint[2*i]]>true[i]),na.rm=TRUE))
417        }
418        est <- rbind(est,Coverage=Coverage)
419    }
420    if (!is.null(names)) {
421         if (length(names)<ncol(est)) {
422            uest <- unique(estimate)
423            names <- names[match(estimate,uest)]
424        }
425        colnames(est) <- names
426
427    }
428    if (unique.names && !is.null(colnames(est))) {
429        colnames(est) <- make.unique(colnames(est))
430    }
431    est[is.nan(est)] <- NA
432
433    return(structure(est,
434                     n=NROW(object),
435                     time=tm,
436                     class=c("summary.sim","matrix")))
437}
438