R version 3.2.4 (2016-03-10) -- "Very Secure Dishes" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## code to reproduce Table 5 > ## THIS FILE: SIMEX simulations for n=5000 > library(splines) > set.seed(4074) > options(width=150) > sessionInfo() R version 3.2.4 (2016-03-10) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] splines stats graphics grDevices utils datasets methods base > objfile <- "~/simex-vfunc/code-files-for-psychometrika-paper-aug2016/code-for-table5-n5000.Robj" > > ################################################ > ## function for deconvolution via NPMLE > ################################################ > > npmle_deconv <- function(W, csem, xmin = -5, xmax = +5, ngrid = 2000, lambda = 0.0005, lltol = 1e-7, psmall = 0.00005, discrete = FALSE, quietly = FALSE){ + ## nonparametric MLE of distribution of error-free variable "X" given + ## error-prone variable "W" where W = X + U and U|X ~ N(0, csem^2(X)) + ## + ## uses algorithm described in Rabe-Hesketh et al (2003). "Correcting for + ## covariate measurement error in logistic regression using nonparametric + ## maximum likelihood estimation." Statistical Modelling 3: 215-232. + ## + ## INPUTS: + ## W: vector of observed scores + ## csem: function of single variable x that returns conditional standard deviation of W|X=x. + ## needs to be able to accept a vector input and return vector output + ## xmin, xmax, ngrid: define grid over which X distribution is calculated + ## NOTE: default (xmin,xmax) are appropriate only for a standardized score; will need to be changed + ## for typical reporting scales for test scores + ## lambda: step in RH algorithm which is on probability scale - this seemed like a good value based on simulations + ## lltol: when we can't find a step that improves the likelihood by more than lltol, we stop + ## psmall: if a point evolves to have mass less than psmall we drop it + ## discrete: TODO for future version + ## quietly: if FALSE, prints iteration-by-iteration output + ## + ## VALUE: + ## list(.history = history of probability distributions, pu = dataframe containing summary info of last distribution) + stopifnot( all(!is.na(W)) && is.function(csem) && (xmin < xmax) && (lambda > 0) && (lltol > 0) ) + if(discrete){ + stop("discrete option to be implemented in future version") + } + if(min(W) < xmin){ + stop("xmin exceeds smallest data value, consider decreasing xmin") + } + if(max(W) > xmax){ + stop("largest data value exceeds xmax; consider increasing xmax") + } + + ## NOTE: we collapse over repeat W values when possible, using "counts" to contribute to likelihood + .W <- sort(unique(W)) + nW <- length(.W) + .counts <- sapply(.W, function(x){ sum(W == x) }) + stopifnot(sum(.counts) == length(W)) + W <- .W + + ## functions to map probabilities to reduced-dimension unconstrained scale, and inverse + theta_to_p <- function(theta){ + if(length(theta) == 0){ + return(1) + } else { + p <- exp(theta) / (1 + sum(exp(theta))) + return(c(p, 1-sum(p))) + } + } + + p_to_theta <- function(p){ + stopifnot( all(p > 0) && (abs((1 - sum(p))) < 1e-9) ) + if(length(p) == 1){ + return(numeric(0)) + } else { + last <- length(p) + return(log(p[1:(last-1)] * (1 + ((1-p[last])/p[last])) )) + } + } + + ## CHECKS: + ## p <- c(0.1, 0.2, 0.4, 0.3); print(ma(p - theta_to_p(p_to_theta(p)))) + ## theta <- c(-3.2, 1.0, -1.0); print(ma(theta - p_to_theta(theta_to_p(theta)))) + ## print(ma(1 - theta_to_p(p_to_theta(1)))) + ## p_to_theta(theta_to_p(numeric(0))) + + ## build (nW x ngrid) matrix of conditional densities if discrete == FALSE, or conditional probabilities if discrete == TRUE + grid <- seq(from=xmin, to=xmax, length=ngrid) + grid.csem <- csem(grid) + if(!discrete){ + fwx <- dnorm( matrix(W,ncol=ngrid,nrow=nW), mean=matrix(grid,ncol=ngrid,nrow=nW,byrow=T), sd=matrix(grid.csem,ncol=ngrid,nrow=nW,byrow=T)) + } else { ## uses pxgu - fwx[i,j] = p(W=W[i] | X = grid[j]) + tmp <- lapply(grid, function(u){ pxgu(u, csem(u), .W)}) + stopifnot(all(sapply(tmp, function(x){ all(x[,1] == .W) }))) + fwx <- do.call("cbind", lapply(tmp, function(x){ x[,2] })) + } + + ## negative log likelihood for a set of probabilities given ".inds" which are + ## indices of "grid" and which are continually updated + ## NOTE: updated to use counts + negll <- function(theta){ + -sum(.counts * log(as.vector(fwx[,.inds] %*% theta_to_p(theta)))) + } + + ## find initial best grid point. NOTE this will often want the X with the + ## biggest CSEM because a single point is inconsistent with most W unless the + ## CSEM is large. the exception is if we pick ugrid to be very large, then it + ## will find interior point. however even if it picks an extreme starting + ## point, that point will be dropped later if it is too far in the tails. + ll <- apply(fwx, 2, function(x){ sum(.counts * log(x)) }) + .inds <- which.max(ll) + .probs <- 1 + .eligible <- setdiff(1:ngrid, .inds) + ll.current <- ll[.inds] + + .history <- list() + .history[[1]] <- data.frame(x = grid[.inds], csem = grid.csem[.inds], p = 1, ll = ll.current, ex = grid[.inds], varx = 0) + + ## now add grid points per the algorithm until there is no improvement + done <- FALSE + while(!done){ + ## evaluate each eligible point with a weight "lambda" + part.static <- as.vector(fwx[,.inds,drop=FALSE] %*% ((1 - lambda)*.probs)) + part.total <- matrix(part.static, ncol=length(.eligible), nrow=nW) + (fwx[,.eligible] * lambda) + ll.candidate <- apply(part.total, 2, function(x){ sum(.counts * log(x)) }) + if(all(ll.candidate - ll.current < lltol)){ + done <- TRUE + } else { + w <- which.max(ll.candidate - ll.current) + .inds <- c(.inds, .eligible[w]) + .eligible <- setdiff(.eligible, .eligible[w]) + + ## set starting value: mass 0.05 at new value, normalized p on existing values + o <- optim(p_to_theta(c(0.95*(.probs/sum(.probs)), 0.05)), negll, method = "BFGS", control = list(maxit=1000, trace=5*(1 - as.numeric(quietly)), REPORT = 1)) + .probs <- theta_to_p(o$par) + ll.current <- -o$value + + ## sometimes we might have picked up an early grid point but now it has virtually no probability, so dump it + w <- which(.probs < psmall) + if(length(w) > 0){ + .probs <- .probs[-w] + .probs <- .probs/sum(.probs) + .inds <- .inds[-w] + } + + ## reorder the grid if need be + o <- order(.inds) + .inds <- .inds[o] + .probs <- .probs[o] + + ## summarize the state + .history[[length(.history)+1]] <- data.frame(x = grid[.inds], + csem = grid.csem[.inds], + p = .probs, + ll = ll.current, + ex = sum(grid[.inds] * .probs), + varx = sum(grid[.inds]^2 * .probs) - (sum(grid[.inds] * .probs))^2) + if(!quietly){ + print(.history[[length(.history)]]) + cat("\n\n\n") + } + } + } + return(list(.history = .history, px = .history[[length(.history)]])) + } > > ################################################ > ## SIMEX simulation > ################################################ > n <- 5000 > nsim <- 500 > B <- 250 > lambda <- seq(from=0, to=2, length=20) > nlambda <- length(lambda) > nparam <- 1 > > bigres <- vector(nsim, mode="list") > > for(simi in 1:nsim){ + #################################################### + ## data generation + #################################################### + + ## generate X bivariate standard normal with correlation 0.4 + d <- data.frame(x1 = rnorm(n)) + + ## generate y according to nonlinear regression (one parameter) + .beta <- 2 + d$egx <- 1 / (1 + exp(-.beta *d$x1)) + d$y <- d$egx + rnorm(n, sd = sqrt( var(d$egx) * (1 - 0.8)/0.8)) ## rsq about 0.8 + + ## generate W as conditionally normal + g1 <- function(x){ + 0.07*(x+1)^2 + } + + d$w1 <- d$x1 + rnorm(n, sd = sqrt(g1(d$x1))) + + ## corresponding A(W) function, closed form in this case + A1 <- function(w){ + (.07/(.07 + 1))*(w+1)^2 + } + + ## get deconvolution estimate of distribution of X. note use rounded W to + ## speed up computations + d$w1round <- round(d$w1, digits=1) + pxhat <- npmle_deconv(d$w1round, csem = function(x){ sqrt(g1(x))}, xmin=min(d$w1round), xmax=max(d$w1round), quietly=TRUE)$px + + ## get pieces we need for various SIMEX + d$egx <- var(d$w1) - unique(pxhat$varx) + + ## conditional mean of X given W (so we can calculate g(E[X|W]) + f <- function(w){ + pxw <- dnorm(w, mean = pxhat$x, sd = sqrt(g1(pxhat$x))) * pxhat$p + pxw <- pxw / sum(pxw) + return( sum(pxhat$x * pxw) ) + } + d$exgivenw <- sapply(d$w1, f) + + ## E[g(X) | W] + f <- function(w){ + pxw <- dnorm(w, mean = pxhat$x, sd = sqrt(g1(pxhat$x))) * pxhat$p + pxw <- pxw / sum(pxw) + return( sum(g1(pxhat$x) * pxw) ) + } + d$egxgivenw <- sapply(d$w1, f) + + ########################################################################### + ## infeasible SIMEX using unobserved true variance for each case + ########################################################################### + rL.gx <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * g1(d$x1))) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.gx[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.gx <- as.data.frame(do.call("rbind",rL.gx)) + + ########################################################################### + ## SIMEX using estimated average ME, ignoring heteroskedasticity + ########################################################################### + rL.avggx <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * d$egx)) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.avggx[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.avggx <- as.data.frame(do.call("rbind",rL.avggx)) + + ########################################################################### + ## SIMEX using naive estimate of average ME variance + ########################################################################### + rL.avggw <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * mean(g1(d$w1)))) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.avggw[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.avggw <- as.data.frame(do.call("rbind",rL.avggw)) + + ########################################################################### + ## naive SIMEX using g1(w) plug-in + ########################################################################### + rL.gw <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * g1(d$w1))) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.gw[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.gw <- as.data.frame(do.call("rbind",rL.gw)) + + ########################################################################### + ## naive SIMEX using g1(E[X|W]) + ########################################################################### + rL.gexw <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * g1(d$exgivenw))) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.gexw[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.gexw <- as.data.frame(do.call("rbind",rL.gexw)) + + ########################################################################### + ## SIMEX using A(W) + ########################################################################### + rL.aw <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * A1(d$w1))) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.aw[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.aw <- as.data.frame(do.call("rbind",rL.aw)) + + ########################################################################### + ## SIMEX using E[g(X)|W] + ########################################################################### + rL.rcal <- vector(nlambda, mode="list") + + for(.l in 1:nlambda){ + tmp <- vector(B, mode="list") + for(.b in 1:B){ + d$w1star <- d$w1 + rnorm(n, sd = sqrt(lambda[.l] * d$egxgivenw)) + tmp[[.b]] <- coef(nls(formula = y ~ 1 / (1 + exp(-b * w1star)), data=d, start = list(b=1))) + } + rL.rcal[[.l]] <- apply(do.call("rbind", tmp), 2, mean) + } + rL.rcal <- as.data.frame(do.call("rbind",rL.rcal)) + + ############################################################################ + ## collect estimates + ############################################################################ + res.ideal <- data.frame(method = "ideal", coef(nls(formula = y ~ 1 / (1 + exp(-b * x1)), data=d, start = list(b=1))), stringsAsFactors=FALSE) + names(res.ideal) <- c("method","b1") + + res.naive <- data.frame(method = "naive", coef(nls(formula = y ~ 1 / (1 + exp(-b * w1)), data=d, start = list(b=1))), stringsAsFactors=FALSE) + names(res.naive) <- c("method","b1") + + res.regcal <- data.frame(method = "regcalib", coef(nls(formula = y ~ 1 / (1 + exp(-b * exgivenw)), data=d, start = list(b=1))), stringsAsFactors=FALSE) + names(res.regcal) <- c("method","b1") + + proj <- list(rL.gx = rL.gx, + rL.avggx = rL.avggx, + rL.avggw = rL.avggw, + rL.gw = rL.gw, + rL.gexw = rL.gexw, + rL.aw = rL.aw, + rL.rcal = rL.rcal) + + ## quadratic projections + res.rest <- do.call("rbind", lapply( proj, function(x){ + e <- rep(0,nparam) + for(i in 1:length(e)){ + ex <- data.frame(lam = lambda, e = x[,i]) + synth <- data.frame(lam = seq(from=-1, to=2, length=1000)) + mod2 <- lm(e ~ lam + I(lam^2), data=ex) + synth$ehat <- predict(mod2, newdata=synth) + e[i] <- synth$ehat[1] + } + e + })) + res.rest <- as.data.frame(res.rest) + names(res.rest) <- c("b1") + res.rest$method <- paste("simex2", sapply(strsplit(rownames(res.rest),"\\."), function(x){ x[2] }), sep="_") + res.rest <- res.rest[,c("method","b1")] + + ## quartic projections + res.rest4 <- do.call("rbind", lapply( proj, function(x){ + e <- rep(0,nparam) + for(i in 1:length(e)){ + ex <- data.frame(lam = lambda, e = x[,i]) + synth <- data.frame(lam = seq(from=-1, to=2, length=1000)) + mod4 <- lm(e ~ lam + I(lam^2) + I(lam^3) + I(lam^4), data=ex) + synth$ehat <- predict(mod4, newdata=synth) + e[i] <- synth$ehat[1] + } + e + })) + res.rest4 <- as.data.frame(res.rest4) + names(res.rest4) <- c("b1") + res.rest4$method <- paste("simex4", sapply(strsplit(rownames(res.rest4),"\\."), function(x){ x[2] }), sep="_") + res.rest4 <- res.rest4[,c("method","b1")] + + ## put the pieces together + res <- rbind(res.ideal, res.naive, res.regcal, res.rest, res.rest4) + res$sim <- simi + bigres[[simi]] <- list(res = res, projections = proj) + print(simi) + } [1] 1 [1] 2 [1] 3 [1] 4 [1] 5 [1] 6 [1] 7 [1] 8 [1] 9 [1] 10 [1] 11 [1] 12 [1] 13 [1] 14 [1] 15 [1] 16 [1] 17 [1] 18 [1] 19 [1] 20 [1] 21 [1] 22 [1] 23 [1] 24 [1] 25 [1] 26 [1] 27 [1] 28 [1] 29 [1] 30 [1] 31 [1] 32 [1] 33 [1] 34 [1] 35 [1] 36 [1] 37 [1] 38 [1] 39 [1] 40 [1] 41 [1] 42 [1] 43 [1] 44 [1] 45 [1] 46 [1] 47 [1] 48 [1] 49 [1] 50 [1] 51 [1] 52 [1] 53 [1] 54 [1] 55 [1] 56 [1] 57 [1] 58 [1] 59 [1] 60 [1] 61 [1] 62 [1] 63 [1] 64 [1] 65 [1] 66 [1] 67 [1] 68 [1] 69 [1] 70 [1] 71 [1] 72 [1] 73 [1] 74 [1] 75 [1] 76 [1] 77 [1] 78 [1] 79 [1] 80 [1] 81 [1] 82 [1] 83 [1] 84 [1] 85 [1] 86 [1] 87 [1] 88 [1] 89 [1] 90 [1] 91 [1] 92 [1] 93 [1] 94 [1] 95 [1] 96 [1] 97 [1] 98 [1] 99 [1] 100 [1] 101 [1] 102 [1] 103 [1] 104 [1] 105 [1] 106 [1] 107 [1] 108 [1] 109 [1] 110 [1] 111 [1] 112 [1] 113 [1] 114 [1] 115 [1] 116 [1] 117 [1] 118 [1] 119 [1] 120 [1] 121 [1] 122 [1] 123 [1] 124 [1] 125 [1] 126 [1] 127 [1] 128 [1] 129 [1] 130 [1] 131 [1] 132 [1] 133 [1] 134 [1] 135 [1] 136 [1] 137 [1] 138 [1] 139 [1] 140 [1] 141 [1] 142 [1] 143 [1] 144 [1] 145 [1] 146 [1] 147 [1] 148 [1] 149 [1] 150 [1] 151 [1] 152 [1] 153 [1] 154 [1] 155 [1] 156 [1] 157 [1] 158 [1] 159 [1] 160 [1] 161 [1] 162 [1] 163 [1] 164 [1] 165 [1] 166 [1] 167 [1] 168 [1] 169 [1] 170 [1] 171 [1] 172 [1] 173 [1] 174 [1] 175 [1] 176 [1] 177 [1] 178 [1] 179 [1] 180 [1] 181 [1] 182 [1] 183 [1] 184 [1] 185 [1] 186 [1] 187 [1] 188 [1] 189 [1] 190 [1] 191 [1] 192 [1] 193 [1] 194 [1] 195 [1] 196 [1] 197 [1] 198 [1] 199 [1] 200 [1] 201 [1] 202 [1] 203 [1] 204 [1] 205 [1] 206 [1] 207 [1] 208 [1] 209 [1] 210 [1] 211 [1] 212 [1] 213 [1] 214 [1] 215 [1] 216 [1] 217 [1] 218 [1] 219 [1] 220 [1] 221 [1] 222 [1] 223 [1] 224 [1] 225 [1] 226 [1] 227 [1] 228 [1] 229 [1] 230 [1] 231 [1] 232 [1] 233 [1] 234 [1] 235 [1] 236 [1] 237 [1] 238 [1] 239 [1] 240 [1] 241 [1] 242 [1] 243 [1] 244 [1] 245 [1] 246 [1] 247 [1] 248 [1] 249 [1] 250 [1] 251 [1] 252 [1] 253 [1] 254 [1] 255 [1] 256 [1] 257 [1] 258 [1] 259 [1] 260 [1] 261 [1] 262 [1] 263 [1] 264 [1] 265 [1] 266 [1] 267 [1] 268 [1] 269 [1] 270 [1] 271 [1] 272 [1] 273 [1] 274 [1] 275 [1] 276 [1] 277 [1] 278 [1] 279 [1] 280 [1] 281 [1] 282 [1] 283 [1] 284 [1] 285 [1] 286 [1] 287 [1] 288 [1] 289 [1] 290 [1] 291 [1] 292 [1] 293 [1] 294 [1] 295 [1] 296 [1] 297 [1] 298 [1] 299 [1] 300 [1] 301 [1] 302 [1] 303 [1] 304 [1] 305 [1] 306 [1] 307 [1] 308 [1] 309 [1] 310 [1] 311 [1] 312 [1] 313 [1] 314 [1] 315 [1] 316 [1] 317 [1] 318 [1] 319 [1] 320 [1] 321 [1] 322 [1] 323 [1] 324 [1] 325 [1] 326 [1] 327 [1] 328 [1] 329 [1] 330 [1] 331 [1] 332 [1] 333 [1] 334 [1] 335 [1] 336 [1] 337 [1] 338 [1] 339 [1] 340 [1] 341 [1] 342 [1] 343 [1] 344 [1] 345 [1] 346 [1] 347 [1] 348 [1] 349 [1] 350 [1] 351 [1] 352 [1] 353 [1] 354 [1] 355 [1] 356 [1] 357 [1] 358 [1] 359 [1] 360 [1] 361 [1] 362 [1] 363 [1] 364 [1] 365 [1] 366 [1] 367 [1] 368 [1] 369 [1] 370 [1] 371 [1] 372 [1] 373 [1] 374 [1] 375 [1] 376 [1] 377 [1] 378 [1] 379 [1] 380 [1] 381 [1] 382 [1] 383 [1] 384 [1] 385 [1] 386 [1] 387 [1] 388 [1] 389 [1] 390 [1] 391 [1] 392 [1] 393 [1] 394 [1] 395 [1] 396 [1] 397 [1] 398 [1] 399 [1] 400 [1] 401 [1] 402 [1] 403 [1] 404 [1] 405 [1] 406 [1] 407 [1] 408 [1] 409 [1] 410 [1] 411 [1] 412 [1] 413 [1] 414 [1] 415 [1] 416 [1] 417 [1] 418 [1] 419 [1] 420 [1] 421 [1] 422 [1] 423 [1] 424 [1] 425 [1] 426 [1] 427 [1] 428 [1] 429 [1] 430 [1] 431 [1] 432 [1] 433 [1] 434 [1] 435 [1] 436 [1] 437 [1] 438 [1] 439 [1] 440 [1] 441 [1] 442 [1] 443 [1] 444 [1] 445 [1] 446 [1] 447 [1] 448 [1] 449 [1] 450 [1] 451 [1] 452 [1] 453 [1] 454 [1] 455 [1] 456 [1] 457 [1] 458 [1] 459 [1] 460 [1] 461 [1] 462 [1] 463 [1] 464 [1] 465 [1] 466 [1] 467 [1] 468 [1] 469 [1] 470 [1] 471 [1] 472 [1] 473 [1] 474 [1] 475 [1] 476 [1] 477 [1] 478 [1] 479 [1] 480 [1] 481 [1] 482 [1] 483 [1] 484 [1] 485 [1] 486 [1] 487 [1] 488 [1] 489 [1] 490 [1] 491 [1] 492 [1] 493 [1] 494 [1] 495 [1] 496 [1] 497 [1] 498 [1] 499 [1] 500 There were 50 or more warnings (use warnings() to see the first 50) > > ## bigres <- do.call("rbind", bigres) > save(bigres, file=objfile, compress=TRUE) > > ## now restrict to summary > bigres <- do.call("rbind", lapply(bigres, function(x){ x$res })) > > ## summaries > for(v in setdiff(names(bigres), c("method","sim"))){ + bigres[,paste(v,"mean",sep="_")] <- ave(bigres[,v], bigres$method, FUN=mean) + bigres[,paste(v,"sd",sep="_")] <- ave(bigres[,v], bigres$method, FUN=sd) + } > > print(sums <- unique(bigres[,c("method",names(bigres)[grep("_",names(bigres))])]), digits=3) method b1_mean b1_sd b ideal 2.00 0.0285 b1 naive 1.75 0.0271 b2 regcalib 1.88 0.0308 rL.gx simex2_gx 2.06 0.0394 rL.avggx simex2_avggx 2.18 0.0402 rL.avggw simex2_avggw 2.21 0.0409 rL.gw simex2_gw 1.95 0.0357 rL.gexw simex2_gexw 2.02 0.0371 rL.aw simex2_aw 1.93 0.0348 rL.rcal simex2_rcal 2.05 0.0376 rL.gx1 simex4_gx 2.06 0.0580 rL.avggx1 simex4_avggx 2.27 0.0589 rL.avggw1 simex4_avggw 2.32 0.0618 rL.gw1 simex4_gw 1.77 0.0402 rL.gexw1 simex4_gexw 1.94 0.0442 rL.aw1 simex4_aw 1.78 0.0410 rL.rcal1 simex4_rcal 1.97 0.0454 > ## true value is > print(.beta) [1] 2 > > > proc.time() user system elapsed 141499.464 17.989 141520.740