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