1## Called from ./lmrob_simulation.Rnw 2## ~~~~~~~~~~~~~~~~~~~~~ 3 4########################################################################### 5## 1. simulation helper functions 6########################################################################### 7 8f.estname <- function(est = 'lmrob') 9 ## Purpose: translate between 'estname' and actual function name, 10 ## defaults to 'lmrob' 11 ## f.lmRob is just a wrapper for lmRob, since there are some 12 ## problems with the weight and weights arguments 13 ## ---------------------------------------------------------------------- 14 ## Arguments: est: name of estimator 15 ## ---------------------------------------------------------------------- 16 ## Author: Manuel Koller, Date: 6 Oct 2009, 13:36 17 switch(est, 18 lm.rbase = 'lmrob', lm.robust = 'f.lmRob', rlm = 'rlm', lm = 'lm', 19 est) 20 21f.errname <- function(err, prefix = 'r') 22 ## Purpose: translate between natural name of distribution and 23 ## R (r,p,q,d)-name 24 ## ---------------------------------------------------------------------- 25 ## Arguments: err: name of distribution 26 ## ---------------------------------------------------------------------- 27 ## Author: Manuel Koller, Date: 6 Oct 2009, 13:36 28 paste(prefix, 29 switch(err,normal="norm", t="t", cauchy="cauchy",cnormal="cnorm", 30 err),sep = '') 31 32f.requires.envir <- function(estname) 33 ## Purpose: returns indicator on whether estname requires envir argument 34 ## ---------------------------------------------------------------------- 35 ## Arguments: estname: name of estimating function 36 ## ---------------------------------------------------------------------- 37 ## Author: Manuel Koller, Date: 7 Oct 2009, 09:34 38 switch(estname, 39 f.lmrob.local = TRUE, 40 FALSE) 41 42f..paste..list <- function(lst) 43 if (length(lst) == 0) return("") else 44 paste(names(lst),lst,sep='=',collapse=', ') 45 46f..split..str <- function(str) { 47 litems <- strsplit(str,', ') 48 lst <- lapply(litems, function(str) strsplit(str,'=')) 49 rlst <- list() 50 for (llst in lst) { 51 lv <- vector() 52 for (litem in llst) lv[litem[1]] <- litem[2] 53 rlst <- c(rlst, list(lv)) 54 } 55 rlst 56} 57 58f.list2str <- function(lst, idx) 59 ## Purpose: convert a list into a string that identifies the 60 ## function and parameter configuration 61 ## ---------------------------------------------------------------------- 62 ## Arguments: lst: list or list of lists 63 ## idx: only take the elements in idx 64 ## ---------------------------------------------------------------------- 65 ## Author: Manuel Koller, Date: 7 Oct 2009, 10:03 66 f..paste..list(if(missing(idx)) unlist(lst) else unlist(lst)[idx]) 67 68f.as.numeric <- function(val) 69{ 70 ## Purpose: convert value to numeric if possible 71 ## ---------------------------------------------------------------------- 72 ## Arguments: vec: value to convert 73 ## ---------------------------------------------------------------------- 74 ## Author: Manuel Koller, Date: 26 Oct 2009, 12:10 75 76 r <- suppressWarnings(as.numeric(val)) 77 if (is.na(r)) { 78 ## is character, try to convert to TRUE and FALSE 79 return(switch(casefold(val), 80 "true" = TRUE, 81 "false" = FALSE, 82 val)) 83 } else return(r) 84} 85f.as.numeric.vectorized <- function(val) sapply(val, f.as.numeric) 86 87f.as.integer <- function(val) 88{ 89 ## Purpose: convert value to numeric if possible 90 ## ---------------------------------------------------------------------- 91 ## Arguments: vec: value to convert 92 ## ---------------------------------------------------------------------- 93 ## Author: Manuel Koller, Date: 26 Oct 2009, 12:10 94 95 r <- suppressWarnings(as.integer(val)) 96 if (is.na(r)) { 97 ## is character, try to convert to TRUE and FALSE 98 return(switch(casefold(val), 99 "true" = TRUE, 100 "false" = FALSE, 101 val)) 102 } else return(r) 103} 104 105f.str2list <- function(str, splitchar = '\\.') 106{ 107 ## Purpose: inverse of f.list2str 108 ## ---------------------------------------------------------------------- 109 ## Arguments: str: string or list of strings produced with f.list2str 110 ## ---------------------------------------------------------------------- 111 ## Author: Manuel Koller, Date: 8 Oct 2009, 14:20 112 113 ## split input string or strings into a list of vectors 114 lst <- f..split..str(as.character(str)) 115 rlst <- list() 116 ## walk list 117 for (lv in lst) { 118 lrlst <- list() 119 ## for each element of the vector 120 for (ln in names(lv)) { 121 ## split 122 lnames <- strsplit(ln, splitchar)[[1]] 123 ## set either directly 124 if (length(lnames) == 1) lrlst[ln] <- f.as.numeric(lv[ln]) 125 ## or, if it contains a dot, as a sublist 126 else { 127 if (is.null(lrlst[[lnames[1]]])) lrlst[[lnames[1]]] <- list() 128 lrlst[[lnames[1]]][paste(lnames[-1],collapse='.')] <- f.as.numeric(lv[ln]) 129 } 130 } 131 rlst <- c(rlst, list(lrlst)) 132 } 133 rlst 134} 135 136f.round.numeric <- function(num, digits = 0) { ## round only numeric values in list 137 idx <- sapply(num, is.numeric) 138 ret <- num 139 ret[idx] <- lapply(num[idx],round,digits=digits) 140 ret 141} 142 143f.errs2str <- function(errs) 144{ 145 ## Purpose: convert list of errors into pretty strings 146 ## ---------------------------------------------------------------------- 147 ## Arguments: errs: estlist element errs 148 ## ---------------------------------------------------------------------- 149 ## Author: Manuel Koller, Date: 8 Oct 2009, 14:51 150 151 rv <- vector() 152 for (lerr in errs) { 153 rv <- c(rv, 154 switch(lerr$err, 155 normal = paste("N(",lerr$args$mean,",", 156 lerr$args$sd,")", sep=""), 157 set =, 158 t = paste("t",lerr$args$df,sep=""), 159 paste(lerr$err,"(",paste(f.round.numeric(lerr$args,2), 160 collapse=","),")",sep=""))) 161 } 162 rv 163} 164 165f.procedures2str <- function(procs) 166{ 167 ## Purpose: convert procedures element in estlist to pretty data.frame 168 ## ---------------------------------------------------------------------- 169 ## Arguments: proc: estlist element procedures 170 ## ---------------------------------------------------------------------- 171 ## Author: Manuel Koller, Date: 8 Oct 2009, 14:57 172 173 rdf <- rep(" ",7) 174 175 for (lproc in procs) { 176 method <- if(is.null(lproc$args$method)) 177 switch(lproc$estname, 178 lm = 'lsq', 179 "SM") else lproc$args$method 180 cov <- switch(lproc$estname, ## lm.robust, rlm, lmrob: set default arguments 181 lm.robust = list(cov = 'Default', 182 cov.corrfact = 'empirical', 183 cov.xwx = TRUE, 184 cov.resid = 'trick', 185 cov.hubercorr = TRUE, 186 cov.dfcorr = 1), 187 rlm = list(cov = 'Default', 188 cov.corrfact = 'empirical', 189 cov.xwx = FALSE, 190 cov.resid = 'final', 191 cov.hubercorr = TRUE, 192 cov.dfcorr = 1), 193 ## lmrob = list(cov = 'f.avar1', ## method .vcov.MM equals f.avar1 194 ## cov.resid = 'final'), 195 lmrob = do.call('lmrob.control', ## get default arguments from lmrob.control 196 lproc$args)[c('cov', 'cov.corrfact', 'cov.xwx', 197 'cov.resid', 'cov.hubercorr', 'cov.dfcorr')], 198 if (is.null(lproc$args)) list(cov = 'Default') else lproc$args) 199 if (is.null(lproc$args$psi)) { 200 psi <- switch(lproc$estname, 201 rlm =, 202 lmrob = 'bisquare', 203 lm.robust = { 204 if (is.null(lproc$args$weight)) { 205 if (is.null(lproc$args$weight2)) 'optimal' 206 else lproc$args$weight2 207 } else lproc$args$weight[2] }, 208 "NA") 209 } else { 210 psi <- lproc$args$psi 211 ## test if tuning.psi is the default one 212 if (!is.null(lproc$args$tuning.psi) && 213 isTRUE(all.equal(lproc$args$tuning.psi, .Mpsi.tuning.default(psi)))) 214 psi <- paste(psi, lproc$args$tuning.psi) 215 } 216 D.type <- switch(lproc$estname, 217 lmrob.u =, 218 lmrob = if (is.null(lproc$args$method) || 219 lproc$args$method %in% c('SM', 'MM')) 'S' else 'D', 220 lmrob.mar = if (is.null(lproc$args$type)) 'qE' else lproc$args$type, 221 rlm = 'rlm', 222 lm.robust = 'rob', 223 lm = 'lm', 224 'NA') 225 rdf <- rbind(rdf,c(lproc$estname, method, f.args2str(lproc$args), 226 cov$cov, f.cov2str(cov), psi, D.type)) 227 } 228 colnames(rdf) <- c("Function", "Method", "Tuning", "Cov", "Cov.Tuning", "Psi", "D.type") 229 if (NROW(rdf) == 2) t(rdf[-1,]) else rdf[-1,] 230} 231 232f.chop <- function(str,l=1) 233 ## Purpose: chop string by l characters 234 ## ---------------------------------------------------------------------- 235 ## Arguments: str: string to chop 236 ## l: number of characters to chop 237 ## ---------------------------------------------------------------------- 238 ## Author: Manuel Koller, Date: 8 Oct 2009, 15:19 239 substr(str,1,nchar(str)-l) 240 241fMpsi2str <- function(psi) 242{ 243 ## Purpose: make pretty M.psi and D.chi, etc. 244 ## ---------------------------------------------------------------------- 245 ## Arguments: M.psi: M.psi argument 246 ## ---------------------------------------------------------------------- 247 ## Author: Manuel Koller, Date: 8 Oct 2009, 15:28 248 249 if (is.null(psi)) psi 250 else if (psi == "tukeyPsi1" || psi == "tukeyChi") "bisquare" 251 else if (grepl("Psi1$",psi)) f.chop(psi,4) 252 else if (grepl("Chi$",psi)) f.chop(psi,3) 253 else psi 254} 255 256f.c.psi2str <- function(c.psi) 257{ 258 ## Purpose: make pretty tuning.psi and D.tuning.chi, etc. 259 ## ---------------------------------------------------------------------- 260 ## Arguments: c.psi: tuning.psi argument 261 ## ---------------------------------------------------------------------- 262 ## Author: Manuel Koller, Date: 8 Oct 2009, 15:34 263 264 if (is.null(c.psi)) return(NULL) 265 266 round(as.numeric(c.psi),2) 267} 268 269f.args2str <- function(args) 270{ 271 ## Purpose: convert args element in procedures element of estlist 272 ## to a pretty string 273 ## ---------------------------------------------------------------------- 274 ## Arguments: args: args element in procedures element of estlist 275 ## ---------------------------------------------------------------------- 276 ## Author: Manuel Koller, Date: 8 Oct 2009, 15:11 277 278 lst <- list() 279 lst$psi <- if (!is.null(args$weight)) args$weight[2] 280 else if (!is.null(args$weight2)) args$weight2 281 else args$psi 282 283 lst$c.psi <- if (!is.null(args$efficiency)) 284 round(f.eff2c.psi(args$efficiency, lst$psi),2) 285 else f.c.psi2str(args$tuning.psi) 286 287 if (!is.null(args$method) && grepl("D",args$method)) { 288 lst$D <- if (!is.null(args$D.type)) args$D.type else NULL 289 lst$tau <- args$tau 290 } 291 292 f..paste..list(lst) 293} 294 295f.cov2str <- function(args) 296{ 297 ## Purpose: convert cov part in args element in procedures element of 298 ## estlist to a pretty string 299 ## ---------------------------------------------------------------------- 300 ## Arguments: args: args element in procedures element of estlist 301 ## ---------------------------------------------------------------------- 302 ## Author: Manuel Koller, Date: 8 Oct 2009, 15:39 303 304 lst <- list() 305 306 if (!is.null(args$cov) && !args$cov %in% c('Default','f.avarwh')) 307 lst$cov <- sub('^f\\.', '', args$cov) 308 else { 309 lst$hc <- args$cov.hubercorr 310 lst$dfc <- args$cov.dfcorr 311 lst$r <- args$cov.resid 312 lst$rtau <- args$cov.corrfact 313 lst$xwx <- args$cov.xwx 314 } 315 ## convert logical to numeric 316 lst <- lapply(lst, function(x) if (is.logical(x)) as.numeric(x) else x) 317 318 f..paste..list(lst) 319} 320 321f.procstr2id <- function(procstrs, fact = TRUE) 322{ 323 ## Purpose: create short identifiers of procstrs 324 ## ---------------------------------------------------------------------- 325 ## Arguments: procstrs: vector of procstrs 326 ## fact: convert to factor or not 327 ## ---------------------------------------------------------------------- 328 ## Author: Manuel Koller, Date: 3 Nov 2009, 08:58 329 330 lst0 <- f.str2list(procstrs) 331 r <- sapply(lst0, function(x) { 332 paste(c(x$estname, 333 if (is.null(x$args$method)) NULL else x$args$method, 334 substr(c(x$args$psi,x$args$weight2, x$args$weight[2]), 1, 3)), 335 collapse = '.') 336 }) 337 if (fact) ru <- unique(r) 338 if (fact) factor(r, levels = ru, labels = ru) else r 339} 340 341f.splitstrs <- function(strs, split = '_', ...) 342{ 343 ## Purpose: split vector of strings by split and convert the list into 344 ## a data.frame with columns type and id 345 ## ---------------------------------------------------------------------- 346 ## Arguments: strs: vector of strings 347 ## split: character vector to use for splitting 348 ## ...: arguments to strsplit, see ?strsplit 349 ## ---------------------------------------------------------------------- 350 ## Author: Manuel Koller, Date: 19 Oct 2009, 08:46 351 352 lstr <- strsplit(strs, split, ...) 353 ldf <- t(as.data.frame(lstr)) 354 rownames(ldf) <- NULL 355 as.data.frame(ldf, stringsAsFactors = FALSE) 356} 357 358f.abind <- function(arr1,arr2, along = ndim) 359{ 360 ## Purpose: like abind, but less powerful 361 ## ---------------------------------------------------------------------- 362 ## Arguments: arr1, arr2: arrays to bind 363 ## along: dimension along to bind to, 364 ## defaults to last dimension 365 ## ---------------------------------------------------------------------- 366 ## Author: Manuel Koller, Date: 20 Oct 2009, 11:33 367 ## if along =! last dimension: permutate array 368 ndim <- length(dim(arr1)) 369 if (along != ndim) { 370 arr1 <- aperm(arr1, perm = c((1:ndim)[-along],along)) 371 arr2 <- aperm(arr2, perm = c((1:ndim)[-along],along)) 372 } 373 ldmn1 <- dimnames(arr1) 374 ldmn2 <- dimnames(arr2) 375 ld1 <- dim(arr1) 376 ld2 <- dim(arr2) 377 if (length(ld1) != length(ld2)) 378 stop('f.abind: Dimensions must be identical') 379 if (!identical(ldmn1[-ndim],ldmn2[-ndim])) 380 stop('f.abind: Dimnames other than in the along dimension must match exactly') 381 if (any(ldmn1[[ndim]] %in% ldmn2[[ndim]])) 382 stop('f.abind: Dimnames in along dimension must be unique') 383 ldmn3 <- ldmn1 384 ldmn3[[ndim]] <- c(ldmn1[[ndim]], ldmn2[[ndim]]) 385 ld3 <- ld1 386 ld3[ndim] <- ld1[ndim] + ld2[ndim] 387 ## build array 388 arr3 <- array(c(arr1, arr2), dim = ld3, dimnames = ldmn3) 389 ## permutate dimensions back 390 if (along != ndim) { 391 lperm <- 1:ndim 392 lperm[along] <- ndim 393 lperm[(along+1):ndim] <- along:(ndim-1) 394 arr3 <- aperm(arr3, perm = lperm) 395 } 396 arr3 397} 398 399f.abind.3 <- function(...) f.abind(..., along = 3) 400 401f.rename.level <- function(factor, from, to) { 402 ## Purpose: rename level in a factor 403 ## ---------------------------------------------------------------------- 404 ## Arguments: factor: factor variable 405 ## from: level to be changed 406 ## to: value 407 ## ---------------------------------------------------------------------- 408 ## Author: Manuel Koller, Date: 18 Aug 2010, 14:45 409 levels(factor)[levels(factor) == from] <- to 410 factor 411} 412 413 414 415########################################################################### 416## 2. main simulation functions 417########################################################################### 418 419f.sim <- function(estlist, 420 .combine = 'f.abind', 421 .combine.2 = 'f.abind.3', ## hack for foreach 422 silent = TRUE) 423{ 424 ## Purpose: perform simulation according to estlist entry ec 425 ## ---------------------------------------------------------------------- 426 ## Arguments: ec: estlist, list consisting of: 427 ## - design: data frame of design 428 ## - nrep: number of repetitions 429 ## - errs: list of error distributions including arguments 430 ## - err: name of error distribution 431 ## - args: list of arguments (to be passed to do.call() 432 ## - procedures: list of parameter configurations and 433 ## procedures to call 434 ## - estname: name of estimation procedure 435 ## - args: arguments that define the call 436 ## silent: silent argument to try 437 ## ---------------------------------------------------------------------- 438 ## Author: Werner Stahel / Manuel Koller, Date: 21 Aug 2008, 07:55 439 440 ## get designs 441 ldd <- estlist$design 442 use.intercept <- if(is.null(estlist$use.intercept)) TRUE 443 else estlist$use.intercept 444 nobs <- NROW(ldd) 445 npar <- NCOL(ldd) + use.intercept 446 nrep <- estlist$nrep 447 nlerrs <- nobs*nrep 448 ## initialize: 449 lestlist <- estlist 450 ## 'evaluate' estlist$procedure list 451 lprocs <- c() 452 for (i in 1:length(estlist$procedures)) { 453 ## generate lprocstr (identification string) 454 lprocs[i] <- estlist[['procedures']][[i]][['lprocstr']] <- 455 f.list2str(estlist[['procedures']][[i]]) 456 } 457 ## find all error distributions 458 lerrs <- unique(sapply(lestlist$errs, f.list2str)) 459 ## walk estlist$output to create output column names vector 460 ## store result into lnames, it is used in f.sim.process 461 lnames <- c() 462 for (i in 1:length(estlist$output)) { 463 llnames <- estlist[['output']][[i]][['lnames']] <- 464 eval(estlist[['output']][[i]][['names']]) 465 lnames <- c(lnames, llnames) 466 } 467 468 ## get different psi functions 469 lpsifuns <- unlist(unique(lt <- sapply(estlist$procedures, function(x) x$args$psi))) 470 ## get entries without psi argument 471 lrest <- sapply(lt, is.null) 472 if (sum(lrest) > 0) lpsifuns <- c(lpsifuns, '__rest__') 473 474 ## Walk error distributions 475 res <- foreach(lerrlst = estlist$errs, .combine = .combine) %:% 476 foreach(lpsifun = lpsifuns, .combine = .combine.2) %dopar% { 477 ## filter for psi functions 478 lidx <- if (lpsifun == '__rest__') lrest else 479 unlist(sapply(estlist$procedures, 480 function(x) !is.null(x$args$psi) && x$args$psi == lpsifun)) 481 cat(f.errs2str(list(lerrlst)), lpsifun, " ") 482 ## get function name and parameters 483 lerrfun <- f.errname(lerrlst$err) 484 lerrpar <- lerrlst$args 485 lerrstr <- f.list2str(lerrlst) 486 487 ## --- initialize array 488 lres <- array(NA, dim=c(nrep, ## data dimension 489 length(lnames), ## output type dimension 490 sum(lidx), ## estimation functions and arguments dimension 491 1), ## error distributions dimension 492 dimnames = list(Data = NULL, 493 Type = lnames, Procstr = lprocs[lidx], Errstr = lerrstr)) 494 ## set seed 495 set.seed(estlist$seed) 496 ## generate errors: seperately for each repetition 497 lerrs <- c(sapply(1:nrep, function(x) do.call(lerrfun, c(n = nobs, lerrpar)))) 498 499 ## if estlist$design has an attribute 'gen' 500 ## then this function gen will generate designs 501 ## and takes arguments: n, p, rep 502 ## and returns the designs in a list 503 if (is.function(attr(ldd, 'gen'))) { 504 ldds <- attr(ldd, 'gen')(nobs, npar - use.intercept, nrep, lerrlst) 505 } 506 507 ## Walk repetitions 508 for (lrep in 1:nrep) { 509 if (lrep%%100 == 0) cat(" ", lrep) 510 lerr <- lerrs[(1:nobs)+(lrep-1)*nobs] 511 if (exists('ldds')) { 512 ldd <- ldds[[lrep]] 513 ## f.sim.reset.envirs() 514 } 515 ## Walk estimator configurations 516 for (lproc in estlist$procedures[lidx]) { 517 ## call estimating procedure 518 lrr <- tryCatch(do.call(f.estname(lproc$estname), 519 c(if(use.intercept) 520 list(lerr ~ . , data = ldd) else 521 list(lerr ~ . - 1, data = ldd), lproc$args)), 522 error=function(e)e) 523 ERR <- inherits(lrr, 'error') 524 if (ERR && !silent) { 525 print(lproc$lprocstr) 526 print(lrr) 527 } 528 if (!silent && !converged(lrr)) { 529 print(lproc$lprocstr) 530 browser() ## <<< 531 } 532 ## check class: if procedure failed: 533 if (ERR) next 534 ## check convergence of estimator 535 if (!converged(lrr)) next 536 ## process output 537 for (lov in estlist$output) { 538 llnames <- lov$lnames 539 ret <- tryCatch(lres[lrep,llnames,lproc$lprocstr,lerrstr] <- eval(lov$fun), 540 error= function(e)e) 541 if (!silent && inherits(ret, 'error')) { 542 cat('Error', dQuote(ret$message), 'in repetition',lrep, 543 '\n for:',llnames,'procstr:',lproc$lprocstr,'\n') 544 browser() ## <<< 545 print(lov$fun) 546 print(try(eval(lov$fun))) 547 } 548 } 549 } 550 } 551 ## print debug information if requested 552 if (!silent) str(lres) 553 lres 554 } 555 ## restore original order of lprocs 556 res <- res[,,match(lprocs, dimnames(res)[[3]]),,drop=FALSE] 557 ## set attributes 558 attr(res, 'estlist') <- lestlist 559 cat("\n") 560 res 561} 562 563########################################################################### 564## build estlist 565########################################################################### 566 567f.combine <- function(..., keep.list = FALSE) { 568 ## Purpose: creates a list of all combinations of elements given as 569 ## arguments, similar to expand.grid. 570 ## Arguments can be named. 571 ## If an argument is a list, then its elements are considered 572 ## as fixed objects that should not be recombined. 573 ## if keep.list = TRUE, these elements are combined 574 ## as a list with argument. 575 ## ---------------------------------------------------------------------- 576 ## Arguments: collection of lists or vectors with argument names 577 ## ---------------------------------------------------------------------- 578 ## Author: Manuel Koller, Date: 7 Oct 2009, 11:13 579 ## convert arguments into a big list 580 args <- list(...) 581 ## if more than two arguments, call recursively 582 if (length(args) > 2) 583 lst <- do.call("f.combine", c(args[-1], list(keep.list=keep.list))) 584 else { 585 ## if just two arguments, create list of second argument 586 ## if this is a list, then there's nothing to do 587 if (!keep.list && is.list(args[[2]])) lst <- args[[2]] 588 ## else convert to a list of one-elements lists with proper name 589 else { 590 lst <- list() 591 for (lelem in args[[2]]) { 592 llst <- list(lelem) 593 if (!is.null(names(args)[2])) names(llst)[1] <- names(args)[2] 594 lst <- c(lst, list(llst)) 595 } 596 } 597 } 598 ## ok, now we can add the first element to all elements of lst 599 lst2 <- list() 600 if (keep.list && is.list(args[[1]])) args[[1]] <- lapply(args[[1]], list) 601 for (lelem in args[[1]]) { 602 for (relem in lst) { 603 llst <- c(lelem, relem) 604 if (nchar(names(llst)[1]) == 0 && nchar(names(args)[1])>0) 605 names(llst)[1] <- names(args)[1] 606 lst2 <- c(lst2, list(llst)) 607 } 608 } 609 lst2 610} 611 612## some fragments to build estlist 613## errors 614.errs.normal.1 <- list(err = 'normal', 615 args = list(mean = 0, sd = 1)) 616.errs.normal.2 <- list(err = 'normal', 617 args = list(mean = 0, sd = 2)) 618.errs.t.13 <- list(err = 't', 619 args = list(df = 13)) 620.errs.t.11 <- list(err = 't', 621 args = list(df = 11)) 622.errs.t.10 <- list(err = 't', 623 args = list(df = 10)) 624.errs.t.9 <- list(err = 't', 625 args = list(df = 9)) 626.errs.t.8 <- list(err = 't', 627 args = list(df = 8)) 628.errs.t.7 <- list(err = 't', 629 args = list(df = 7)) 630.errs.t.5 <- list(err = 't', 631 args = list(df = 5)) 632.errs.t.3 <- list(err = 't', 633 args = list(df = 3)) 634.errs.t.1 <- list(err = 't', 635 args = list(df = 1)) 636 637## skewed t distribution 638.errs.skt.Inf.2 <- list(err = 'cskt', 639 args = list(df = Inf, gamma = 2)) 640.errs.skt.5.2 <- list(err = 'cskt', 641 args = list(df = 5, gamma = 2)) 642## log normal distribution 643.errs.lnrm <- list(err = 'lnorm', 644 args = list(meanlog = 0, sdlog = 0.6936944)) 645## laplace distribution 646.errs.laplace <- list(err = 'laplace', 647 args = list(location = 0, scale = 1/sqrt(2))) 648 649## contaminated normal 650.errs.cnorm..1.0.10 <- list(err = 'cnorm', 651 args = list(epsilon = 0.1, meanc = 0, sdc = sqrt(10))) 652 653.errs.cnorm..1.4.1 <- list(err = 'cnorm', 654 args = list(epsilon = 0.1, meanc = 4, sdc = 1)) 655 656.errs.test <- list(.errs.normal.1 657 ,.errs.t.5 658 ,.errs.t.3 659 ,.errs.t.1 660 ) 661 662## arguments 663.args.final <- f.combine(psi = c('optimal', 'bisquare', 'lqq', 'hampel'), 664 seed = 0, 665 max.it = 500, 666 k.max = 2000, 667 c(list(list(method = 'MM', cov = '.vcov.avar1')), 668 list(list(method = 'MM', cov = '.vcov.w', 669 start = 'lrr')), 670 f.combine(method = c('SMD', 'SMDM'), 671 cov = '.vcov.w', 672 start = 'lrr'))) 673 674## use fixInNamespace("lmrob.fit", "robustbase") 675## insert: 676## N = { 677## tmp <- lmrob..M..fit(x = x/init$tau, y = y/init$tau, obj = 678## init) 679## tmp$qr <- NULL 680## tmp 681## }, 682 683## .args.final <- f.combine(psi = c('optimal', 'bisquare', 'ggw', 'lqq'), 684## seed = 0, 685## max.it = 500, 686## k.max = 2000, 687## c(list(list(method = "SMDM", cov = '.vcov.w')), 688## list(list(method = "SMDN", cov = '.vcov.w', 689## start = 'lrr')))) 690 691 692## standard for lmRob 693.args.bisquare.lmRob.0 <- list(## initial.alg = 'random', 694 efficiency = 0.95 695 ,weight = c('bisquare', 'bisquare'), 696 trace = FALSE 697 ) 698 699.args.optimal.lmRob.0 <- list(## initial.alg = 'random', 700 efficiency = 0.95 701 ,weight = c('optimal', 'optimal'), 702 trace = FALSE) 703 704.procedures.final <- c(list(list(estname = 'lm')), 705 f.combine(estname = 'lmrob.u', args = .args.final, 706 keep.list = TRUE), 707 f.combine(estname = 'lmrob.mar', 708 args = f.combine(psi = 'bisquare', 709 seed = 0, max.it = 500, k.max = 2000, 710 cov = '.vcov.w', type = c('qT', 'qE')), 711 keep.list = TRUE), 712 f.combine(estname = 'lm.robust', 713 args = list(.args.bisquare.lmRob.0, 714 .args.optimal.lmRob.0), keep.list = TRUE)) 715 716## output 717.output.sigma <- list(sigma = list( 718 names = quote("sigma"), 719 fun = quote(sigma(lrr)))) 720.output.beta <- list(beta = list( 721 names = quote(paste('beta',1:npar,sep='_')), 722 fun = quote(coef(lrr)))) 723.output.se <- list(se = list( 724 names = quote(paste('se',1:npar,sep='_')), 725 fun = quote(sqrt(diag(covariance.matrix(lrr)))))) 726.output.sumw <- list(sumw = list( 727 names = quote("sumw"), 728 fun = quote(sum(robustness.weights(lrr))))) 729.output.nnz <- list(nnz = list( 730 names = quote("nnz"), 731 fun = quote(sum(robustness.weights(lrr) < 1e-3)))) 732 733########################################################################### 734## simulation results processing functions 735########################################################################### 736 737## use apply to aggregate data 738## use matplot(t(result)) to plot aggregated data 739 740f.apply <- function(res, items = dimnames(res)[[2]], 741 FUN, ..., swap = FALSE) 742{ 743 ## Purpose: similar to apply, return data not as matrix, but 744 ## as data.frame 745 ## ---------------------------------------------------------------------- 746 ## Arguments: res: simulation results array 747 ## items: items to use in apply 748 ## FUN: function to apply 749 ## ...: additional arguments to FUN 750 ## swap: if TRUE: swap first two columns 751 ## ---------------------------------------------------------------------- 752 ## Author: Manuel Koller, Date: 8 Oct 2009, 13:39 753 754 ## aggregate data 755 lz <- apply(res[,items,,,drop=FALSE], 2:4, FUN, ...) 756 ## if return object has four dimensions (multidim output of FUN) 757 ## rotate first three dimensions 758 if (length(dim(lz)) == 4 && swap) aperm(lz, perm=c(2,1,3,4)) else lz 759} 760 761f.dimnames2df <- function(arr, dm = dimnames(arr), 762 page = TRUE, err.on.same.page = TRUE, 763 value.col = ndim - 2, 764 procstr.col = ndim - 1, 765 errstr.col = ndim, 766 procstr.id = TRUE, 767 split = '_') 768{ 769 ## Purpose: create data frame from dimnames: 770 ## len_1 .. len_100, cpr_1 .. cpr_100 771 ## will yield a data frame with column id from 1 .. 100 772 ## column type with cpr and len and columns procstr and errstr 773 ## It is assumed, that the max number (100) is the same for all 774 ## output value types 775 ## ---------------------------------------------------------------------- 776 ## Arguments: arr: 3 or more dim array (optional) 777 ## dm: dimnames to be used 778 ## page: add a column page to simplify plots 779 ## err.on.same.page: whether all errs should be on the same 780 ## page 781 ## value.col: index of value column (set to NULL for none) 782 ## the values in this column are split name_id 783 ## and put into two columns in the data frame 784 ## procstr.col: index of procedure column 785 ## (both: or NULL for not to be converted) 786 ## errstr.col: index of error string column 787 ## procstr.id: create procstr id 788 ## ---------------------------------------------------------------------- 789 ## Author: Manuel Koller, Date: 19 Oct 2009, 08:41 790 791 if (!is.list(dm)) stop('f.dimnames2df: dm must be a list') 792 ## remove 'NULL' dimensions 793 dm <- dm[!sapply(dm,is.null)] 794 ndim <- length(dm) 795 if (ndim == 0) stop('f.dimnames2df: dimnames all null') 796 ldims <- sapply(dm, length) 797 ## split and convert types into data.frame 798 if (!is.null(value.col)) { 799 ldf <- f.splitstrs(dm[[value.col]], split = split) 800 lid <- NCOL(ldf) == 2 801 if (lid) lids <- unique(as.numeric(ldf[,2])) ## convert ids into numeric 802 ## we do not need to repeat over different types of values, only ids 803 ldims[value.col] <- ldims[value.col] / length(unique(ldf[,1])) 804 } 805 ## merge into one large data.frame: for each distribution 806 rdf <- list() 807 for (ld in 1:ndim) { 808 lname <- if (is.null(lname <- names(dm)[ld])) length(rdf)+1 else lname 809 ltimes <- if (ld == ndim) 1 else prod(ldims[(ld+1):ndim]) 810 leach <- if (ld == 1) 1 else prod(ldims[1:(ld-1)]) 811 if (!is.null(value.col) && ld == value.col) { 812 if (lid) rdf[[paste(lname,'Id')]] <- 813 rep(lids,times=ltimes,each=leach) ## value ids 814 ## no else: the values will be added in the a2df procedures 815 } else if (!is.null(procstr.col) && ld == procstr.col) { 816 ## convert procstrs to data.frame with pretty names 817 lprdf <- data.frame(f.procedures2str(f.str2list(dm[[ld]])), 818 Procstr = factor(dm[[ld]], levels = dm[[ld]], 819 labels = dm[[ld]])) 820 if (procstr.id) lprdf$PId <- f.procstr2id(dm[[ld]]) 821 ## repeat 822 lprdf <- if (ltimes == 1 && leach == 1) 823 lprdf else apply(lprdf,2,rep,times=ltimes,each=leach) 824 lprdf <- as.data.frame(lprdf, stringsAsFactors=FALSE) 825 ## convert all into nice factors (with the original ordering) 826 for (lk in colnames(lprdf)) { 827 luniq <- unique(lprdf[[lk]]) 828 lprdf[[lk]] <- factor(lprdf[[lk]], levels = luniq, labels = luniq) 829 } 830 rdf <- c(rdf, lprdf) 831 } else if (!is.null(errstr.col) && ld == errstr.col) { 832 ## convert errstrs to data.frame with pretty names 833 ledf <- f.errs2str(f.str2list(dm[[ld]])) 834 ## repeat and convert to factor with correct ordering 835 rdf[[lname]] <- factor(rep(dm[[ld]],times=ltimes,each=leach), 836 levels = dm[[ld]], labels = dm[[ld]]) 837 rdf[['Error']] <- factor(rep(ledf,times=ltimes,each=leach), 838 levels = ledf, labels = ledf) 839 } else { 840 ## no conversion necessary 841 rdf[[lname]] <- rep(dm[[ld]],times=ltimes,each=leach) 842 } 843 } 844 ## add page argument 845 if (page && !is.null(procstr.col)) { 846 ltpf <- if (!is.null(errstr.col) && !err.on.same.page) 847 interaction(rdf[['Procstr']],rdf[['Error']]) 848 else interaction(rdf[['Procstr']]) 849 rdf[['Page']] <- as.numeric(factor(ltpf, unique(ltpf))) 850 } 851 rdf <- as.data.frame(rdf) 852 if (!is.null(value.col)) 853 attr(rdf, 'Types') <- unique(ldf[,1]) 854 rdf 855} 856 857f.a2df.2 <- function(arr, dm = dimnames(arr), err.on.same.page = FALSE, ...) 858{ 859 ## Purpose: convert arr to data.frame 860 ## uses f.dimnames2df and adds a column to contain the values 861 ## if ndim == 4 and dimnames NULL: assumes first dimension is 862 ## data dimension which is ignored by f.dimnames2df 863 ## add counter 864 ## ---------------------------------------------------------------------- 865 ## Arguments: arr: array to convert 866 ## ---------------------------------------------------------------------- 867 ## Author: Manuel Koller, Date: 23 Oct 2009, 12:29 868 869 ## ndim == 2 ?? 870 ndim <- length(dim(arr)) 871 ## if ndim == 4: check if dimnames of dim 1 are NULL 872 if (ndim == 4 && is.null(dm[[1]])) 873 dm[[1]] <- 1:dim(arr)[1] 874 rdf <- f.dimnames2df(dm=dm, ...) 875 ## just add values for all 'Types', possibly including Type.ID 876 if (ndim > 2) 877 for (lvt in attr(rdf, 'Types')) { 878 llvt <- if (is.null(rdf$Type.Id)) lvt else paste(lvt,unique(rdf$Type.Id),sep='_') 879 rdf[[lvt]] <- as.vector(switch(ndim, 880 stop('wrong number of dimensions'), ## 1 881 arr, ## 2 882 arr[llvt,,], ## 3 883 arr[,llvt,,])) ## 4 884 } 885 else 886 rdf$values <- as.vector(arr) 887 rdf 888} 889 890 891f.dimnames2pc.df <- function(arr, dm = dimnames(arr), 892 npcs = NCOL(estlist$design.predict), ...) 893{ 894 ## Purpose: create data frame to be used in plotting of pc components 895 ## calls f.dimnames2df and adds an additional column for 896 ## identifying the principal components 897 ## ---------------------------------------------------------------------- 898 ## Arguments: arr, dm: see f.dimnames.df 899 ## npcs: number of principal components 900 ## ---------------------------------------------------------------------- 901 ## Author: Manuel Koller, Date: 23 Oct 2009, 11:51 902 903 if (missing(npcs) && !is.null(attr(estlist$design.predict, 'npcs'))) 904 npcs <- attr(estlist$design.predict, 'npcs') 905 906 ## convert into data.frame 907 rdf <- f.dimnames2df(dm = dm, ...) 908 ## calculate number of points per principal component 909 npts <- (length(unique(rdf$Type.Id)) - 1) / npcs 910 ## add new column pc 911 rdf$PC <- 1 912 if (npcs > 1) 913 for (li in 2:npcs) { 914 lids <- (1:npts + npts*(li-1) + 1) 915 rdf$PC[rdf$Type.Id %in% lids] <- li ## fixme: center is not repeated 916 } 917 rdf$PC <- factor(rdf$PC, levels = 1:npcs, labels = paste('PC',1:npcs,sep=' ')) 918 rdf 919} 920 921f.a2pc.df <- function(arr, ...) 922{ 923 ## Purpose: convert arr to data.frame 924 ## uses f.dimnames2pc.df and adds a column to contain the values 925 ## ---------------------------------------------------------------------- 926 ## Arguments: arr: array to convert 927 ## ---------------------------------------------------------------------- 928 ## Author: Manuel Koller, Date: 23 Oct 2009, 12:29 929 930 ## convert dimnames 931 rdf <- f.dimnames2pc.df(arr, err.on.same.page = FALSE,...) 932 ## add values 933 for (lvt in attr(rdf, 'Types')) 934 rdf[[lvt]] <- as.vector(arr[paste(lvt,unique(rdf$Type.Id),sep='_'),,]) 935 ## repeat values: only PC_1 has center value, add it for other PCs 936 ## build index 937 idx <- 1:NROW(rdf) 938 rpc <- as.character(rdf$PC) 939 for (lerr in levels(rdf$Error)) { 940 for (lprc in levels(rdf$Procstr)) { 941 for (lpc in levels(rdf$PC)) { 942 if (lpc == 'PC 1') next 943 ## get first entry of this PC 944 lmin <- min(which(rdf$Error == lerr & rdf$Procstr == lprc & rdf$PC == lpc)) 945 ## where is this in idx? 946 lwm <- min(which(lmin == idx)) 947 ## get first entry of PC_1 948 lmin1 <- min(which(rdf$Error == lerr & rdf$Procstr == lprc & rdf$PC == 'PC 1')) 949 ## update idx 950 idx <- c(idx[1:(lwm-1)], lmin1, idx[lwm:length(idx)]) 951 ## update PC column of result 952 rpc <- c(rpc[1:(lwm-1)], lpc, rpc[lwm:length(rpc)]) 953 } 954 } 955 } 956 ## repeat centers 957 rdf <- rdf[idx,] 958 ## update PC column 959 rdf$PC <- factor(rpc) 960 ## return 961 rdf 962} 963 964f.calculate <- function(expr,arr,dimname = as.character(expr)) 965{ 966 ## Purpose: calculate formula and return as conformable array 967 ## ---------------------------------------------------------------------- 968 ## Arguments: expr: expression to calculate (string is also ok) 969 ## arr: array (from f.sim) 970 ## dimname: name of the calculated dimension 971 ## ---------------------------------------------------------------------- 972 ## Author: Manuel Koller, Date: 9 Oct 2009, 10:15 973 974 if (!is.expression(expr)) expr <- as.expression(expr) 975 976 lnams <- dimnames(arr)[[2]] 977 lst <- list() 978 for (lnam in lnams) 979 expr <- gsub(paste(lnam,'\\b',sep=''), 980 paste("arr[,",lnam,",,,drop=FALSE]",sep='"'), expr) 981 982 r <- eval(parse(text = expr)) 983 dimnames(r)[[2]] <- dimname 984 r 985 986 ## maybe use abind to merge the two arrays? 987} 988 989f.calculate.many <- function(expr, arr, dimname = dims, dims) 990{ 991 ## Purpose: calculate formula and abind into array 992 ## supply expr as string with # symbols to be replaced 993 ## dimname can also contain # symbols 994 ## ---------------------------------------------------------------------- 995 ## Arguments: same as f.calculate and 996 ## dims: vector of items to replace # with 997 ## ---------------------------------------------------------------------- 998 ## Author: Manuel Koller, Date: 14 Oct 2009, 10:11 999 1000 for (i in 1:length(dims)) { 1001 lexpr <- gsub("#",dims[i],expr) 1002 ldimname <- 1003 if (length(dimname) > 1) dimname[i] else gsub("#",dims[i],dimname) 1004 if (i == 1) 1005 rarr <- f.calculate(lexpr,arr,ldimname) 1006 else 1007 rarr <- abind(rarr, f.calculate(lexpr,arr,ldimname), along=2) 1008 } 1009 1010 rarr 1011} 1012 1013f.errs <- function(estlist, err, rep, gen = NULL, nobs, npar) 1014{ 1015 ## Purpose: generate and return errors of specified repetition 1016 ## or, if missing, all errors as a matrix 1017 ## ---------------------------------------------------------------------- 1018 ## Arguments: estlist: estlist 1019 ## err: error distribution (estlist$errs[1] for example) 1020 ## rep: desired repetition (optional) 1021 ## gen: function to generate designs (optional) 1022 ## nobs: nr. rows, npap: nr. predictors (both optional) 1023 ## --------------------------------------------------------------------- 1024 ## Author: Manuel Koller, Date: 13 Oct 2009, 11:21 1025 1026 nobs <- NROW(estlist$design) 1027 nrep <- estlist$nrep 1028 nlerrs <- nobs*nrep 1029 npred <- NROW(estlist$design.predict) 1030 1031 ## get function name and parameters 1032 lerrfun <- f.errname(err$err) 1033 lerrpar <- err$args 1034 lerrstr <- f.list2str(err) 1035 ## set seed 1036 set.seed(estlist$seed) 1037 ## generate errors: seperately for each repetition 1038 lerrs <- c(sapply(1:nrep, function(x) do.call(lerrfun, c(n = nobs, lerrpar)))) 1039 ## lerrs <- do.call(lerrfun, c(n = nlerrs, lerrpar)) 1040 ## to get to the same seed state as f.sim(.default) 1041 ## generate also the additional errors 1042 ## calculate additional number of errors 1043 for (i in 1:length(estlist$output)) { 1044 if (!is.null(estlist[['output']][[i]][['nlerrs']])) 1045 nlerrs <- nlerrs + eval(estlist[['output']][[i]][['nlerrs']]) 1046 } 1047 if (length(lerrs) < nlerrs) 1048 nowhere <- do.call(lerrfun, c(n = nlerrs - length(lerrs), lerrpar)) 1049 ## generate designs 1050 if (!is.null(gen) && is.function(gen)) { 1051 ldds <- gen(nobs, npar, nrep, err) 1052 } 1053 ## return errors 1054 ret <- if (!missing(rep)) lerrs[1:nobs+(rep-1)*nobs] else matrix(lerrs, nobs) 1055 if (exists('ldds')) attr(ret, 'designs') <- if (!missing(rep)) ldds[[i]] else ldds 1056 ret 1057} 1058 1059f.selection <- function(procstrs = dimnames(r.test)[[3]], 1060 what = c('estname', 'args.method', 'args.psi', 'args.tuning.psi', 1061 'args.type', 'args.weight2', 'args.efficiency'), 1062 restr = '') 1063{ 1064 ## Purpose: get selection of results: first one of the specified estimates 1065 ## ---------------------------------------------------------------------- 1066 ## Arguments: procstrs: what is the selection 1067 ## what: named vector to use in grep 1068 ## restr: do not select estimators with procstr 1069 ## that match this regexp parameters 1070 ## ---------------------------------------------------------------------- 1071 ## Author: Manuel Koller, Date: 2 Nov 2009, 09:06 1072 1073 ## match restrictions 1074 lrestr <- -(lall <- 1:length(procstrs)) ## no restrictions 1075 if (!missing(restr)) { 1076 lrestr <- grep(restr, procstrs) 1077 if (length(lrestr) == 0) lrestr <- -lall 1078 procstrs <- procstrs[-lrestr] 1079 } 1080 ## procstr2list, but do not split into sublists 1081 lproclst <- f.str2list(procstrs, splitchar='_____') 1082 ## helper function: select only items that occur what 1083 tfun <- function(x) x[what] 1084 lproclst <- lapply(lproclst, tfun) 1085 ## convert back to string 1086 lprocstr <- sapply(lproclst, f.list2str) 1087 ## get all unique combinations and the first positions 1088 lidx <- match(unique(lprocstr), lprocstr) 1089 r <- procstrs[lidx] 1090 attr(r, 'idx') <- lall[-lrestr][lidx] 1091 r 1092} 1093 1094f.get.current.dimnames <- function(i,dn,margin) 1095{ 1096 ## Purpose: get current dimnames in the margins of array 1097 ## we're applying on 1098 ## ---------------------------------------------------------------------- 1099 ## Arguments: i: counter 1100 ## dn: dimnames 1101 ## margin: margin argument to apply 1102 ## ---------------------------------------------------------------------- 1103 ## Author: Manuel Koller, Date: 16 Apr 2010, 10:44 1104 1105 ## pos <- integer(0) 1106 lcdn <- character(0) 1107 1108 for (lm in margin) { 1109 ## get length of current margin 1110 llen <- length(dn[[lm]]) 1111 ## i modulo llen gives the current position in this dimension 1112 lpos <- (if (i > 0) i-1 else 0) %% llen + 1 1113 ## update pos 1114 ## pos <- c(pos, lpos) 1115 ## update lcdn 1116 lcdn <- c(lcdn, dn[[lm]][lpos]) 1117 ## update i: subtract lpos and divide by llen 1118 i <- (i - lpos) / llen + 1 1119 } 1120 lcdn 1121} 1122 1123f.n <- Vectorize(function(design) 1124{ 1125 ## Purpose: get n obs of design 1126 ## ---------------------------------------------------------------------- 1127 ## Arguments: design: design to get n of 1128 ## ---------------------------------------------------------------------- 1129 ## Author: Manuel Koller, Date: 19 Apr 2010, 11:19 1130 1131 NROW(get(design)) 1132}) 1133 1134f.p <- Vectorize(function(design) 1135{ 1136 ## Purpose: get p par of design 1137 ## ---------------------------------------------------------------------- 1138 ## Arguments: design: design to get p of 1139 ## ---------------------------------------------------------------------- 1140 ## Author: Manuel Koller, Date: 19 Apr 2010, 11:19 1141 1142 NCOL(get(design)) + 1 1143}) 1144 1145f.which.min <- function(x, nr = 1) { 1146 ## Purpose: get the indices of the minimal nr of observations 1147 ## ---------------------------------------------------------------------- 1148 ## Arguments: x: vector of values 1149 ## nr: number of indices to return 1150 ## ---------------------------------------------------------------------- 1151 ## Author: Manuel Koller, Date: 4 May 2010, 12:18 1152 match(sort(x)[1:nr], x) 1153} 1154 1155f.which.max <- function(x, nr = 1) f.which.min(-x, nr) 1156 1157## f.get.scale <- function(procstr, proclst = f.str2list(procstr)) 1158## { 1159## ## Purpose: get scale estimate used for procstrs 1160## ## ---------------------------------------------------------------------- 1161## ## Arguments: procstr: procstrs (dimnames(r.test)[[3]]) as output by 1162## ## f.list2str() 1163## ## proclst: list of procedures, as in estlist$procedures 1164## ## ---------------------------------------------------------------------- 1165## ## Author: Manuel Koller, Date: 9 Sep 2010, 13:52 1166 1167## ret <- list() 1168 1169## for (lproc in proclst) { 1170## if (lproc$estname == 'lm') { 1171## ## least squares 1172## ret <- c(ret, list(list(fun='f.lsq'))) 1173## } else { 1174## ## default (S-scale): 1175## fun <- 'lmrob.mscale' 1176## lidx <- names(lproc$args)[na.omit(match(c('psi', 'tuning.chi', 'seed'), 1177## names(lproc$args)))] 1178 1179## if (!is.null(lproc$args$method) && 1180## substr(lproc$args$method,1,3) == 'SMD') { 1181## ## D-scale 1182## fun <- 'lmrob.dscale' 1183## lidx <- names(lproc$args)[na.omit(match(c('psi', 'tuning.psi'), 1184## names(lproc$args)))] 1185## } else if (lproc$estname == 'lmrob.mar' ### continue here 1186## ret <- c(ret, list(list(fun=fun, args=lproc$args[lidx]))) 1187## } 1188 1189## }) 1190 1191 1192 1193########################################################################### 1194## functions related to prediction 1195########################################################################### 1196 1197f.prediction.points <- function(design, type = c('pc', 'grid'), 1198 length.out = 4*NCOL(design), f = 0.5, 1199 direction = +1, max.pc = 5) 1200{ 1201 ## Purpose: generate prediction points for design 1202 ## generate four points along the second principal component 1203 ## in the center, 2 intermediate distances and long distance 1204 ## (from the center) 1205 ## ---------------------------------------------------------------------- 1206 ## Arguments: design: design matrix 1207 ## type: type of prediction points: grid / principal components 1208 ## length.out: approximate number of prediction points 1209 ## f: extend range by f (like extendrange()) 1210 ## direction: +1 or -1: which direction to go from the center 1211 ## max.pc: maximum number of principal components to use 1212 ## ---------------------------------------------------------------------- 1213 ## Author: Manuel Koller, Date: 9 Oct 2009, 16:48 1214 1215 ## match type argument 1216 type = match.arg(type) 1217 ## get ranges 1218 lrange <- apply(design, 2, range) 1219 ## extend range by f 1220 lrange <- data.frame(apply(lrange, 2, extendrange, f = f)) 1221 1222 switch(type, 1223 pc = { 1224 ## calculate robust covariance matrix 1225 rob <- covMcd(design) 1226 ## and use it to calculate the principal components 1227 rpc <- princomp(covmat = rob$cov) 1228 ## get corner with maximum distance from rob$center 1229 lidx <- apply(abs(lrange - rob$center),2,which.max) 1230 lcr <- diag(as.matrix(lrange[lidx,])) 1231 ## create grid points: 1232 rdf <- rob$center 1233 ## for each principal component 1234 for (id in 1:min(NCOL(rpc$loadings),max.pc)) { 1235 ## calculate factor to reach each boundary 1236 lfct <- (lcr - rob$center) / rpc$loadings[,id] 1237 ## calculate distances to boundaries and take the minimal one 1238 lmin <- which.min(sapply(lfct, function(x) sum((rpc$loadings[,id] * x)^2))) 1239 ## create sequence of multiplicands 1240 lmult <- seq(0,lfct[lmin],length.out=length.out/NCOL(rpc$loadings)) 1241 rdf <- rbind(rdf, rep(rob$center,each=length(lmult)-1) + 1242 direction*lmult[-1] %*% t(rpc$loadings[,id])) 1243 } 1244 }, 1245 grid = { 1246 ## generate sequences for every dimension 1247 lval <- as.data.frame(apply(lrange,2,f.seq, 1248 length.out = round(length.out^(1/NCOL(design))) )) 1249 ## return if 1 dimension, otherwise create all combinations 1250 rdf <- if (NCOL(design) > 1) 1251 t(as.data.frame(do.call('f.combine', lval))) else lval 1252 1253 }) 1254 rdf <- as.data.frame(rdf) 1255 rownames(rdf) <- NULL 1256 colnames(rdf) <- colnames(design) 1257 if (type == 'pc') attr(rdf, 'npcs') <- id 1258 rdf 1259} 1260 1261## ## plot with 1262## require(rgl) 1263## plot3d(design) 1264## points3d(f.prediction.points(design), col = 2) 1265 1266## d.data <- data.frame(y = rnorm(10), x = 1:10) 1267## pred <- f.prediction.points(d.data[,-1,drop=FALSE]) 1268## obj <- f.lmrob.local(y ~ x, d.data) 1269## f.predict(obj, pred, interval = 'prediction') 1270## as.vector(t(cbind(rnorm(4), f.predict(obj, pred, interval = 'prediction')))) 1271 1272## estlist for prediction: 1273## start with .output.test 1274## we only need sigma 1275.output.prediction <- c(.output.sigma,.output.beta,.output.se,.output.sumw,.output.nnz) 1276.output.prediction$predict <- 1277 list(names = quote({ 1278 npred <- NROW(estlist$design.predict) 1279 paste(c('fit', 'lwr', 'upr', 'se.fit', 'cpr'), 1280 rep(1:npred,each = 5), sep = '_')}), 1281 fun = quote({ 1282 lpr <- f.predict(lrr, estlist$design.predict, interval = 'prediction', 1283 se.fit = TRUE) ##, df = 16) 1284 lpr <- cbind(lpr$fit, lpr$se.fit) 1285 lqf <- f.errname(lerrlst$err, 'p') 1286 lcpr <- do.call(lqf, c(list(lpr[,'upr']), lerrpar)) - 1287 do.call(lqf, c(list(lpr[,'lwr']), lerrpar)) 1288 as.vector(t(cbind(lpr,lcpr)))})) 1289 1290.estlist.prediction <- list(design = dd, 1291 nrep = 200, 1292 errs = .errs.test, 1293 seed = 0, 1294 procedures = .procedures.final, 1295 design.predict = f.prediction.points(dd), 1296 output = .output.prediction, 1297 use.intercept = TRUE) 1298 1299## predict confidence intervals instead of prediction intervals 1300.estlist.confint <- .estlist.prediction 1301.estlist.confint$output$predict$fun <- 1302 parse(text=gsub('prediction', 'confidence', deparse(.output.prediction$predict$fun))) 1303 1304########################################################################### 1305## Generate designs - function 1306########################################################################### 1307 1308f.gen <- function(n, p, rep, err) { 1309 ## get function name and parameters 1310 lerrfun <- f.errname(err$err) 1311 lerrpar <- err$args 1312 ## generate random predictors 1313 ret <- lapply(1:rep, function(...) 1314 data.frame(matrix(do.call(lerrfun, c(n = n*p, lerrpar)), n, p))) 1315 attr(ret[[1]], 'gen') <- f.gen 1316 ret 1317} 1318 1319.output.sigmaE <- list(sigmaE = list( 1320 names = quote("sigmaE"), 1321 fun = quote({ 1322 ## estimate scale using current scale estimate. 1323 ## this amounts to recalculating the estimate 1324 ## with just an intercept 1325 llargs <- lproc$args 1326 llestname <- lproc$estname 1327 ## save time and just calculate S-estimate and no covariance matrix 1328 if (grepl('^lmrob', llestname)) { 1329 llestname <- 'lmrob' 1330 llargs$cov <- 'none' 1331 llargs$envir <- NULL ## drop envir argument 1332 if (llargs$method %in% c('MM', 'SM')) llargs$method <- 'S' 1333 if (grepl('M$', llargs$method)) 1334 llargs$method <- f.chop(llargs$method) 1335 } else if (lproc$estname == 'lm.robust') { 1336 llargs$estim <- 'Initial' 1337 } 1338 llrr <- tryCatch(do.call(f.estname(lproc$estname), 1339 c(list(lerr ~ 1), llargs)), 1340 error = function(e)e) 1341 ## check class: if procedure failed: class == 'try-error' 1342 if (inherits(llrr, 'error')) NA 1343 ## check convergence of estimator 1344 else if (lproc$estname != 'lm.robust' && !converged(llrr)) NA 1345 else sigma(llrr) 1346 }))) 1347 1348