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=1000 > library(splines) > set.seed(4079) > 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-n1000-v2.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 <- 1000 > 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 Warning messages: 1: In log(as.vector(fwx[, .inds] %*% theta_to_p(theta))) : NaNs produced 2: In log(as.vector(fwx[, .inds] %*% theta_to_p(theta))) : NaNs produced 3: In log(as.vector(fwx[, .inds] %*% theta_to_p(theta))) : NaNs produced > > ## 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.0578 b1 naive 1.75 0.0592 b2 regcalib 1.88 0.0682 rL.gx simex2_gx 2.06 0.0881 rL.avggx simex2_avggx 2.18 0.0886 rL.avggw simex2_avggw 2.21 0.0900 rL.gw simex2_gw 1.95 0.0789 rL.gexw simex2_gexw 2.02 0.0822 rL.aw simex2_aw 1.93 0.0775 rL.rcal simex2_rcal 2.04 0.0837 rL.gx1 simex4_gx 2.05 0.1284 rL.avggx1 simex4_avggx 2.27 0.1321 rL.avggw1 simex4_avggw 2.31 0.1341 rL.gw1 simex4_gw 1.76 0.0931 rL.gexw1 simex4_gexw 1.93 0.1032 rL.aw1 simex4_aw 1.77 0.0891 rL.rcal1 simex4_rcal 1.97 0.1016 > ## true value is > print(.beta) [1] 2 > > > proc.time() user system elapsed 45134.858 3.938 45118.708