R help - Apr 2011 - How to Extract Information from SIMEX Output

Below is a SIMEX object that was generated with the "simex" function
from the
"simex" package applied to a logistic regression fit. From this
mountain of
information I would like to extract all of the values summarized in this
line:

  .. ..$ variance.jackknife: num [1:5, 1:4] 1.684 1.144 0.85 0.624 0.519 ...

Can someone suggest how to go about doing this? I can extract the upper
level results like fit.simex$coefficients but I have had no success getting
at the lower levels.

Tom

> str(fit.simex)
List of 24
 $ coefficients             : Named num [1:2] -17.1 3
  ..- attr(*, "names")= chr [1:2] "(Intercept)"
"x"
 $ SIMEX.estimates          : num [1:6, 1:3] -1 0 0.5 1 1.5 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "lambda" "(Intercept)" "x"
 $ lambda                   : num [1:5] 0 0.5 1 1.5 2
 $ model                    :List of 31
  ..$ coefficients     : Named num [1:2] -13.27 2.32
  .. ..- attr(*, "names")= chr [1:2] "(Intercept)"
"x"
  ..$ residuals        : Named num [1:1615] -1.12 -1.42 -1.23 -1.07 -1.44
...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ fitted.values    : Named num [1:1615] 0.1032 0.2952 0.1847 0.0656
0.3062 ...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ effects          : Named num [1:1615] 19.552 -9.275 -0.473 -0.283
-0.641 ...
  .. ..- attr(*, "names")= chr [1:1615] "(Intercept)"
"x" "" "" ...
  ..$ R                : num [1:2, 1:2] -15.6 0 -81 -4
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  ..$ rank             : int 2
  ..$ qr               :List of 5
  .. ..$ qr   : num [1:1615, 1:2] -15.6232 0.0292 0.0248 0.0159 0.0295 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:1615] "1" "2" "3"
"4" ...
  .. .. .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..$ rank : int 2
  .. ..$ qraux: num [1:2] 1.02 1.02
  .. ..$ pivot: int [1:2] 1 2
  .. ..$ tol  : num 1e-11
  .. ..- attr(*, "class")= chr "qr"
  ..$ family           :List of 12
  .. ..$ family    : chr "binomial"
  .. ..$ link      : chr "logit"
  .. ..$ linkfun   :function (mu)  
  .. ..$ linkinv   :function (eta)  
  .. ..$ variance  :function (mu)  
  .. ..$ dev.resids:function (y, mu, wt)  
  .. ..$ aic       :function (y, n, mu, wt, dev)  
  .. ..$ mu.eta    :function (eta)  
  .. ..$ initialize:  expression({     if (NCOL(y) == 1) {         if
(is.factor(y))              y <- y != levels(y)[1L]         n <-
rep.int(1,
nobs)         y[weights == 0] <- 0         if (any(y < 0 | y > 1))
stop("y values must be 0 <= y <= 1")         mustart <-
(weights * y +
0.5)/(weights + 1)         m <- weights * y         if (any(abs(m -
round(m)) > 0.001))              warning("non-integer #successes in a
binomial glm!")     }     else if (NCOL(y) == 2) {         if (any(abs(y -
round(y)) > 0.001))              warning("non-integer counts in a
binomial
glm!")         n <- y[, 1] + y[, 2]         y <- ifelse(n == 0, 0,
y[, 1]/n)
weights <- weights * n         mustart <- (n * y + 0.5)/(n + 1)     }    
else stop("for the binomial family, y must be a vector of 0 and
1's\n",
"or a 2 column matrix where col 1 is no. successes and col 2 is no.
failures") })
  .. ..$ validmu   :function (mu)  
  .. ..$ valideta  :function (eta)  
  .. ..$ simulate  :function (object, nsim)  
  .. ..- attr(*, "class")= chr "family"
  ..$ linear.predictors: Named num [1:1615] -2.162 -0.87 -1.485 -2.656
-0.818 ...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ deviance         : num 1521
  ..$ aic              : num 1525
  ..$ null.deviance    : num 1622
  ..$ iter             : int 4
  ..$ weights          : Named num [1:1615] 0.0926 0.208 0.1506 0.0613
0.2125 ...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ prior.weights    : Named num [1:1615] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ df.residual      : int 1613
  ..$ df.null          : int 1614
  ..$ y                : Named num [1:1615] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
  ..$ converged        : logi TRUE
  ..$ boundary         : logi FALSE
  ..$ model            :'data.frame':	1615 obs. of  2 variables:
  .. ..$ y: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1
...
  .. ..$ x: num [1:1615] 4.79 5.35 5.09 4.58 5.37 ...
  .. ..- attr(*, "terms")=Classes 'terms', 'formula'
length 3 y ~ x
  .. .. .. ..- attr(*, "variables")= language list(y, x)
  .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2] "y" "x"
  .. .. .. .. .. ..$ : chr "x"
  .. .. .. ..- attr(*, "term.labels")= chr "x"
  .. .. .. ..- attr(*, "order")= int 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment:
R_GlobalEnv>
  .. .. .. ..- attr(*, "predvars")= language list(y, x)
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2]
"factor" "numeric"
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "y"
"x"
  ..$ x                : num [1:1615, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1615] "1" "2" "3"
"4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ call             : language glm(formula = y ~ x, family = binomial, x
= TRUE, y = TRUE)
  ..$ formula          :Class 'formula' length 3 y ~ x
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  ..$ terms            :Classes 'terms', 'formula' length 3 y ~
x
  .. .. ..- attr(*, "variables")= language list(y, x)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "y" "x"
  .. .. .. .. ..$ : chr "x"
  .. .. ..- attr(*, "term.labels")= chr "x"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(y, x)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "factor"
"numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
  ..$ data             :<environment: R_GlobalEnv> 
  ..$ offset           : NULL
  ..$ control          :List of 3
  .. ..$ epsilon: num 1e-08
  .. ..$ maxit  : num 25
  .. ..$ trace  : logi FALSE
  ..$ method           : chr "glm.fit"
  ..$ contrasts        : NULL
  ..$ xlevels          : Named list()
  ..- attr(*, "class")= chr [1:2] "glm" "lm"
 $ mc.matrix                :List of 1
  ..$ y: num [1:2, 1:2] 0.95 0.05 0.03 0.97
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "0" "1"
  .. .. ..$ : chr [1:2] "0" "1"
 $ B                        : num 800
 $ extrapolation            :List of 12
  ..$ coefficients : num [1:3, 1:2] -13.258 3.299 -0.501 2.315 -0.593 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  ..$ residuals    : num [1:5, 1:2] -0.01292 0.02805 -0.00663 -0.01922
0.01071 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "1" "2" "3" "4"
...
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  ..$ effects      : num [1:5, 1:2] 23.9505 3.63099 -0.46906 -0.00164
0.03846 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "(Intercept)" "lambda"
"I(lambda^2)" "" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  ..$ rank         : int 3
  ..$ fitted.values: num [1:5, 1:2] -13.26 -11.73 -10.46 -9.44 -8.67 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "1" "2" "3" "4"
...
  .. .. ..$ : chr [1:2] "(Intercept)" "x"
  ..$ assign       : int [1:3] 0 1 2
  ..$ qr           :List of 5
  .. ..$ qr   : num [1:5, 1:3] -2.236 0.447 0.447 0.447 0.447 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
  .. .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
  .. .. ..- attr(*, "assign")= int [1:3] 0 1 2
  .. ..$ qraux: num [1:3] 1.45 1.12 1.78
  .. ..$ pivot: int [1:3] 1 2 3
  .. ..$ tol  : num 1e-07
  .. ..$ rank : int 3
  .. ..- attr(*, "class")= chr "qr"
  ..$ df.residual  : int 2
  ..$ xlevels      : Named list()
  ..$ call         : language lm(formula = estimates ~ lambda + I(lambda^2))
  ..$ terms        :Classes 'terms', 'formula' length 3
estimates ~ lambda +
I(lambda^2)
  .. .. ..- attr(*, "variables")= language list(estimates, lambda,
I(lambda^2))
  .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:3] "estimates" "lambda"
"I(lambda^2)"
  .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
  .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
  .. .. ..- attr(*, "order")= int [1:2] 1 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: 0x0bd89124> 
  .. .. ..- attr(*, "predvars")= language list(estimates, lambda,
I(lambda^2))
  .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.2" "numeric"
"numeric"
  .. .. .. ..- attr(*, "names")= chr [1:3] "estimates"
"lambda"
"I(lambda^2)"
  ..$ model        :'data.frame':	5 obs. of  3 variables:
  .. ..$ estimates  : num [1:5, 1:2] -13.27 -11.71 -10.47 -9.46 -8.65 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : NULL
  .. .. .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..$ lambda     : num [1:5] 0 0.5 1 1.5 2
  .. ..$ I(lambda^2):Class 'AsIs'  num [1:5] 0 0.25 1 2.25 4
  .. ..- attr(*, "terms")=Classes 'terms', 'formula'
length 3 estimates ~
lambda + I(lambda^2)
  .. .. .. ..- attr(*, "variables")= language list(estimates, lambda,
I(lambda^2))
  .. .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:3] "estimates" "lambda"
"I(lambda^2)"
  .. .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
  .. .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
  .. .. .. ..- attr(*, "order")= int [1:2] 1 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: 0x0bd89124>
  .. .. .. ..- attr(*, "predvars")= language list(estimates, lambda,
I(lambda^2))
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.2" "numeric"
"numeric"
  .. .. .. .. ..- attr(*, "names")= chr [1:3] "estimates"
"lambda"
"I(lambda^2)"
  ..- attr(*, "class")= chr [1:2] "mlm" "lm"
 $ fitting.method           : chr "quad"
 $ SIMEXvariable            : chr "y"
 $ call                     : language mcsimex(model = fit.naive,
SIMEXvariable = "y", mc.matrix = P,      lambda = c(0.5, 1, 1.5, 2), B
800, fitting.method = "quadratic",  ...
 $ theta                    :List of 2
  ..$ (Intercept):'data.frame':	800 obs. of  4 variables:
  .. ..$ X1: num [1:800] -11.9 -12.4 -11.1 -11.4 -11 ...
  .. ..$ X2: num [1:800] -11.8 -10.05 -9.78 -10.97 -11.12 ...
  .. ..$ X3: num [1:800] -9 -8.71 -9.88 -8.73 -8.65 ...
  .. ..$ X4: num [1:800] -8.49 -8.7 -9.33 -8.23 -8.89 ...
  ..$ x          :'data.frame':	800 obs. of  4 variables:
  .. ..$ X1: num [1:800] 2.07 2.16 1.93 1.97 1.9 ...
  .. ..$ X2: num [1:800] 2.07 1.73 1.69 1.9 1.94 ...
  .. ..$ X3: num [1:800] 1.56 1.49 1.72 1.49 1.48 ...
  .. ..$ X4: num [1:800] 1.47 1.51 1.63 1.42 1.54 ...
 $ fitted.values            : Named num [1:1615] 0.0645 0.2685 0.1421 0.0351
0.2821 ...
  ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
 $ residuals                : Named num [1:1615] -0.0645 -0.2685 -0.1421
-0.0351 -0.2821 ...
  ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
 $ extrapolation.variance   :List of 12
  ..$ coefficients : num [1:3, 1:4] 1.668 -1.105 0.268 -0.321 0.209 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
  .. .. ..$ : NULL
  ..$ residuals    : num [1:5, 1:4] 0.0159 -0.0385 0.0199 0.012 -0.0093 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "1" "2" "3" "4"
...
  .. .. ..$ : NULL
  ..$ effects      : num [1:5, 1:4] -2.1556 -0.9011 0.2503 -0.0154 -0.0461
...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "(Intercept)" "lambda"
"I(lambda^2)" "" ...
  .. .. ..$ : NULL
  ..$ rank         : int 3
  ..$ fitted.values: num [1:5, 1:4] 1.668 1.182 0.83 0.612 0.528 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:5] "1" "2" "3" "4"
...
  .. .. ..$ : NULL
  ..$ assign       : int [1:3] 0 1 2
  ..$ qr           :List of 5
  .. ..$ qr   : num [1:5, 1:3] -2.236 0.447 0.447 0.447 0.447 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
  .. .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
  .. .. ..- attr(*, "assign")= int [1:3] 0 1 2
  .. ..$ qraux: num [1:3] 1.45 1.12 1.78
  .. ..$ pivot: int [1:3] 1 2 3
  .. ..$ tol  : num 1e-07
  .. ..$ rank : int 3
  .. ..- attr(*, "class")= chr "qr"
  ..$ df.residual  : int 2
  ..$ xlevels      : Named list()
  ..$ call         : language lm(formula = variance.jackknife ~ lambda +
I(lambda^2))
  ..$ terms        :Classes 'terms', 'formula' length 3
variance.jackknife ~
lambda + I(lambda^2)
  .. .. ..- attr(*, "variables")= language list(variance.jackknife,
lambda,
I(lambda^2))
  .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:3] "variance.jackknife" "lambda"
"I(lambda^2)"
  .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
  .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
  .. .. ..- attr(*, "order")= int [1:2] 1 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: 0x0bd89124> 
  .. .. ..- attr(*, "predvars")= language list(variance.jackknife,
lambda,
I(lambda^2))
  .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.4" "numeric"
"numeric"
  .. .. .. ..- attr(*, "names")= chr [1:3]
"variance.jackknife" "lambda"
"I(lambda^2)"
  ..$ model        :'data.frame':	5 obs. of  3 variables:
  .. ..$ variance.jackknife: num [1:5, 1:4] 1.684 1.144 0.85 0.624 0.519 ...
  .. ..$ lambda            : num [1:5] 0 0.5 1 1.5 2
  .. ..$ I(lambda^2)       :Class 'AsIs'  num [1:5] 0 0.25 1 2.25 4
  .. ..- attr(*, "terms")=Classes 'terms', 'formula'
length 3
variance.jackknife ~ lambda + I(lambda^2)
  .. .. .. ..- attr(*, "variables")= language list(variance.jackknife,
lambda, I(lambda^2))
  .. .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:3] "variance.jackknife"
"lambda" "I(lambda^2)"
  .. .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
  .. .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
  .. .. .. ..- attr(*, "order")= int [1:2] 1 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: 0x0bd89124>
  .. .. .. ..- attr(*, "predvars")= language list(variance.jackknife,
lambda, I(lambda^2))
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.4" "numeric"
"numeric"
  .. .. .. .. ..- attr(*, "names")= chr [1:3]
"variance.jackknife" "lambda"
"I(lambda^2)"
  ..- attr(*, "class")= chr [1:2] "mlm" "lm"
 $ variance.jackknife       : num [1:2, 1:2] 3.04 -0.58 -0.58 0.111
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..$ : chr [1:2] "(Intercept)" "x"
 $ variance.jackknife.lambda: num [1:6, 1:5] -1 0 0.5 1 1.5 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:6] "1" "" "" "" ...
  .. ..$ : NULL
 $ PSI                      : num [1:1615, 1:10] -0.1032 -0.2952 -0.1847
-0.0656 -0.3062 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:1615] "1" "2" "3" "4"
...
  .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
 $ c11                      : num [1:10, 1:10] 0.151 0.78 0.144 0.748 0.139
...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
  .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
 $ a11                      : num [1:10, 1:10] -0.151 -0.784 0 0 0 ...
 $ sigma                    : num [1:10, 1:10] 2979 -572 2482 -478 2132 ...
 $ sigma.gamma              : num [1:6, 1:6] 2961 -1030 192 -569 195 ...
 $ g                        : num [1:6, 1:2] 1 -1 1 0 0 0 0 0 0 1 ...
 $ s                        : num [1:6, 1:10] -1 0 0 0 0 0 0 0 0 -1 ...
 $ variance.asymptotic      : num [1:2, 1:2] 3.697 -0.704 -0.704 0.134
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "(Intercept)" "x"
  .. ..$ : chr [1:2] "(Intercept)" "x"
 - attr(*, "class")= chr "mcsimex"

--
View this message in context:
http://r.789695.n4.nabble.com/How-to-Extract-Information-from-SIMEX-Output-tp3459082p3459082.html
Sent from the R help mailing list archive at Nabble.com.

On 2011-04-18 17:52, Tom La Bone wrote:
> Below is a SIMEX object that was generated with the "simex"
function from the
> "simex" package applied to a logistic regression fit. From this
mountain of
> information I would like to extract all of the values summarized in this
> line:
>
>    .. ..$ variance.jackknife: num [1:5, 1:4] 1.684 1.144 0.85 0.624 0.519
...
>
> Can someone suggest how to go about doing this? I can extract the upper
> level results like fit.simex$coefficients but I have had no success getting
> at the lower levels.
>
> Tom
This should do it:

  fit.simex$extrapolation.variance$model$variance.jackknife

or

  z <- fit.simex[['extrapolation.variance']]
  z[['model']][['variance.jackknife']]

Peter Ehlers
>
>
>> str(fit.simex)
> List of 24
>   $ coefficients             : Named num [1:2] -17.1 3
>    ..- attr(*, "names")= chr [1:2] "(Intercept)"
"x"
>   $ SIMEX.estimates          : num [1:6, 1:3] -1 0 0.5 1 1.5 ...
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : NULL
>    .. ..$ : chr [1:3] "lambda" "(Intercept)"
"x"
>   $ lambda                   : num [1:5] 0 0.5 1 1.5 2
>   $ model                    :List of 31
>    ..$ coefficients     : Named num [1:2] -13.27 2.32
>    .. ..- attr(*, "names")= chr [1:2] "(Intercept)"
"x"
>    ..$ residuals        : Named num [1:1615] -1.12 -1.42 -1.23 -1.07 -1.44
> ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ fitted.values    : Named num [1:1615] 0.1032 0.2952 0.1847 0.0656
> 0.3062 ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ effects          : Named num [1:1615] 19.552 -9.275 -0.473 -0.283
> -0.641 ...
>    .. ..- attr(*, "names")= chr [1:1615] "(Intercept)"
"x" "" "" ...
>    ..$ R                : num [1:2, 1:2] -15.6 0 -81 -4
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    ..$ rank             : int 2
>    ..$ qr               :List of 5
>    .. ..$ qr   : num [1:1615, 1:2] -15.6232 0.0292 0.0248 0.0159 0.0295 ...
>    .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. ..$ : chr [1:1615] "1" "2" "3"
"4" ...
>    .. .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. ..$ rank : int 2
>    .. ..$ qraux: num [1:2] 1.02 1.02
>    .. ..$ pivot: int [1:2] 1 2
>    .. ..$ tol  : num 1e-11
>    .. ..- attr(*, "class")= chr "qr"
>    ..$ family           :List of 12
>    .. ..$ family    : chr "binomial"
>    .. ..$ link      : chr "logit"
>    .. ..$ linkfun   :function (mu)
>    .. ..$ linkinv   :function (eta)
>    .. ..$ variance  :function (mu)
>    .. ..$ dev.resids:function (y, mu, wt)
>    .. ..$ aic       :function (y, n, mu, wt, dev)
>    .. ..$ mu.eta    :function (eta)
>    .. ..$ initialize:  expression({     if (NCOL(y) == 1) {         if
> (is.factor(y))              y<- y != levels(y)[1L]         n<-
rep.int(1,
> nobs)         y[weights == 0]<- 0         if (any(y<  0 | y>  1))
> stop("y values must be 0<= y<= 1")         mustart<-
(weights * y +
> 0.5)/(weights + 1)         m<- weights * y         if (any(abs(m -
> round(m))>  0.001))              warning("non-integer #successes in
a
> binomial glm!")     }     else if (NCOL(y) == 2) {         if
(any(abs(y -
> round(y))>  0.001))              warning("non-integer counts in a
binomial
> glm!")         n<- y[, 1] + y[, 2]         y<- ifelse(n == 0, 0,
y[, 1]/n)
> weights<- weights * n         mustart<- (n * y + 0.5)/(n + 1)     }
> else stop("for the binomial family, y must be a vector of 0 and
1's\n",
> "or a 2 column matrix where col 1 is no. successes and col 2 is no.
> failures") })
>    .. ..$ validmu   :function (mu)
>    .. ..$ valideta  :function (eta)
>    .. ..$ simulate  :function (object, nsim)
>    .. ..- attr(*, "class")= chr "family"
>    ..$ linear.predictors: Named num [1:1615] -2.162 -0.87 -1.485 -2.656
> -0.818 ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ deviance         : num 1521
>    ..$ aic              : num 1525
>    ..$ null.deviance    : num 1622
>    ..$ iter             : int 4
>    ..$ weights          : Named num [1:1615] 0.0926 0.208 0.1506 0.0613
> 0.2125 ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ prior.weights    : Named num [1:1615] 1 1 1 1 1 1 1 1 1 1 ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ df.residual      : int 1613
>    ..$ df.null          : int 1614
>    ..$ y                : Named num [1:1615] 0 0 0 0 0 0 0 0 0 0 ...
>    .. ..- attr(*, "names")= chr [1:1615] "1"
"2" "3" "4" ...
>    ..$ converged        : logi TRUE
>    ..$ boundary         : logi FALSE
>    ..$ model            :'data.frame':   1615 obs. of  2 variables:
>    .. ..$ y: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1
1 1 1 ...
>    .. ..$ x: num [1:1615] 4.79 5.35 5.09 4.58 5.37 ...
>    .. ..- attr(*, "terms")=Classes 'terms',
'formula' length 3 y ~ x
>    .. .. .. ..- attr(*, "variables")= language list(y, x)
>    .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
>    .. .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. .. ..$ : chr [1:2] "y" "x"
>    .. .. .. .. .. ..$ : chr "x"
>    .. .. .. ..- attr(*, "term.labels")= chr "x"
>    .. .. .. ..- attr(*, "order")= int 1
>    .. .. .. ..- attr(*, "intercept")= int 1
>    .. .. .. ..- attr(*, "response")= int 1
>    .. .. .. ..- attr(*, ".Environment")=<environment:
R_GlobalEnv>
>    .. .. .. ..- attr(*, "predvars")= language list(y, x)
>    .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2]
"factor" "numeric"
>    .. .. .. .. ..- attr(*, "names")= chr [1:2] "y"
"x"
>    ..$ x                : num [1:1615, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:1615] "1" "2" "3"
"4" ...
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. ..- attr(*, "assign")= int [1:2] 0 1
>    ..$ call             : language glm(formula = y ~ x, family = binomial,
x
> = TRUE, y = TRUE)
>    ..$ formula          :Class 'formula' length 3 y ~ x
>    .. .. ..- attr(*, ".Environment")=<environment:
R_GlobalEnv>
>    ..$ terms            :Classes 'terms', 'formula' length
3 y ~ x
>    .. .. ..- attr(*, "variables")= language list(y, x)
>    .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
>    .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. ..$ : chr [1:2] "y" "x"
>    .. .. .. .. ..$ : chr "x"
>    .. .. ..- attr(*, "term.labels")= chr "x"
>    .. .. ..- attr(*, "order")= int 1
>    .. .. ..- attr(*, "intercept")= int 1
>    .. .. ..- attr(*, "response")= int 1
>    .. .. ..- attr(*, ".Environment")=<environment:
R_GlobalEnv>
>    .. .. ..- attr(*, "predvars")= language list(y, x)
>    .. .. ..- attr(*, "dataClasses")= Named chr [1:2]
"factor" "numeric"
>    .. .. .. ..- attr(*, "names")= chr [1:2] "y"
"x"
>    ..$ data             :<environment: R_GlobalEnv>
>    ..$ offset           : NULL
>    ..$ control          :List of 3
>    .. ..$ epsilon: num 1e-08
>    .. ..$ maxit  : num 25
>    .. ..$ trace  : logi FALSE
>    ..$ method           : chr "glm.fit"
>    ..$ contrasts        : NULL
>    ..$ xlevels          : Named list()
>    ..- attr(*, "class")= chr [1:2] "glm" "lm"
>   $ mc.matrix                :List of 1
>    ..$ y: num [1:2, 1:2] 0.95 0.05 0.03 0.97
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:2] "0" "1"
>    .. .. ..$ : chr [1:2] "0" "1"
>   $ B                        : num 800
>   $ extrapolation            :List of 12
>    ..$ coefficients : num [1:3, 1:2] -13.258 3.299 -0.501 2.315 -0.593 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    ..$ residuals    : num [1:5, 1:2] -0.01292 0.02805 -0.00663 -0.01922
> 0.01071 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    ..$ effects      : num [1:5, 1:2] 23.9505 3.63099 -0.46906 -0.00164
> 0.03846 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "(Intercept)" "lambda"
"I(lambda^2)" "" ...
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    ..$ rank         : int 3
>    ..$ fitted.values: num [1:5, 1:2] -13.26 -11.73 -10.46 -9.44 -8.67 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    ..$ assign       : int [1:3] 0 1 2
>    ..$ qr           :List of 5
>    .. ..$ qr   : num [1:5, 1:3] -2.236 0.447 0.447 0.447 0.447 ...
>    .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
>    .. .. ..- attr(*, "assign")= int [1:3] 0 1 2
>    .. ..$ qraux: num [1:3] 1.45 1.12 1.78
>    .. ..$ pivot: int [1:3] 1 2 3
>    .. ..$ tol  : num 1e-07
>    .. ..$ rank : int 3
>    .. ..- attr(*, "class")= chr "qr"
>    ..$ df.residual  : int 2
>    ..$ xlevels      : Named list()
>    ..$ call         : language lm(formula = estimates ~ lambda +
I(lambda^2))
>    ..$ terms        :Classes 'terms', 'formula' length 3
estimates ~ lambda +
> I(lambda^2)
>    .. .. ..- attr(*, "variables")= language list(estimates,
lambda,
> I(lambda^2))
>    .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
>    .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. ..$ : chr [1:3] "estimates" "lambda"
"I(lambda^2)"
>    .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
>    .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
>    .. .. ..- attr(*, "order")= int [1:2] 1 1
>    .. .. ..- attr(*, "intercept")= int 1
>    .. .. ..- attr(*, "response")= int 1
>    .. .. ..- attr(*, ".Environment")=<environment:
0x0bd89124>
>    .. .. ..- attr(*, "predvars")= language list(estimates,
lambda,
> I(lambda^2))
>    .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.2" "numeric"
> "numeric"
>    .. .. .. ..- attr(*, "names")= chr [1:3] "estimates"
"lambda"
> "I(lambda^2)"
>    ..$ model        :'data.frame':       5 obs. of  3 variables:
>    .. ..$ estimates  : num [1:5, 1:2] -13.27 -11.71 -10.47 -9.46 -8.65 ...
>    .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. ..$ : NULL
>    .. .. .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. ..$ lambda     : num [1:5] 0 0.5 1 1.5 2
>    .. ..$ I(lambda^2):Class 'AsIs'  num [1:5] 0 0.25 1 2.25 4
>    .. ..- attr(*, "terms")=Classes 'terms',
'formula' length 3 estimates ~
> lambda + I(lambda^2)
>    .. .. .. ..- attr(*, "variables")= language list(estimates,
lambda,
> I(lambda^2))
>    .. .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
>    .. .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. .. ..$ : chr [1:3] "estimates" "lambda"
"I(lambda^2)"
>    .. .. .. .. .. ..$ : chr [1:2] "lambda"
"I(lambda^2)"
>    .. .. .. ..- attr(*, "term.labels")= chr [1:2]
"lambda" "I(lambda^2)"
>    .. .. .. ..- attr(*, "order")= int [1:2] 1 1
>    .. .. .. ..- attr(*, "intercept")= int 1
>    .. .. .. ..- attr(*, "response")= int 1
>    .. .. .. ..- attr(*, ".Environment")=<environment:
0x0bd89124>
>    .. .. .. ..- attr(*, "predvars")= language list(estimates,
lambda,
> I(lambda^2))
>    .. .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.2" "numeric"
> "numeric"
>    .. .. .. .. ..- attr(*, "names")= chr [1:3]
"estimates" "lambda"
> "I(lambda^2)"
>    ..- attr(*, "class")= chr [1:2] "mlm" "lm"
>   $ fitting.method           : chr "quad"
>   $ SIMEXvariable            : chr "y"
>   $ call                     : language mcsimex(model = fit.naive,
> SIMEXvariable = "y", mc.matrix = P,      lambda = c(0.5, 1, 1.5,
2), B > 800, fitting.method = "quadratic",  ...
>   $ theta                    :List of 2
>    ..$ (Intercept):'data.frame': 800 obs. of  4 variables:
>    .. ..$ X1: num [1:800] -11.9 -12.4 -11.1 -11.4 -11 ...
>    .. ..$ X2: num [1:800] -11.8 -10.05 -9.78 -10.97 -11.12 ...
>    .. ..$ X3: num [1:800] -9 -8.71 -9.88 -8.73 -8.65 ...
>    .. ..$ X4: num [1:800] -8.49 -8.7 -9.33 -8.23 -8.89 ...
>    ..$ x          :'data.frame': 800 obs. of  4 variables:
>    .. ..$ X1: num [1:800] 2.07 2.16 1.93 1.97 1.9 ...
>    .. ..$ X2: num [1:800] 2.07 1.73 1.69 1.9 1.94 ...
>    .. ..$ X3: num [1:800] 1.56 1.49 1.72 1.49 1.48 ...
>    .. ..$ X4: num [1:800] 1.47 1.51 1.63 1.42 1.54 ...
>   $ fitted.values            : Named num [1:1615] 0.0645 0.2685 0.1421
0.0351
> 0.2821 ...
>    ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
>   $ residuals                : Named num [1:1615] -0.0645 -0.2685 -0.1421
> -0.0351 -0.2821 ...
>    ..- attr(*, "names")= chr [1:1615] "1" "2"
"3" "4" ...
>   $ extrapolation.variance   :List of 12
>    ..$ coefficients : num [1:3, 1:4] 1.668 -1.105 0.268 -0.321 0.209 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
>    .. .. ..$ : NULL
>    ..$ residuals    : num [1:5, 1:4] 0.0159 -0.0385 0.0199 0.012 -0.0093
...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. ..$ : NULL
>    ..$ effects      : num [1:5, 1:4] -2.1556 -0.9011 0.2503 -0.0154 -0.0461
> ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "(Intercept)" "lambda"
"I(lambda^2)" "" ...
>    .. .. ..$ : NULL
>    ..$ rank         : int 3
>    ..$ fitted.values: num [1:5, 1:4] 1.668 1.182 0.83 0.612 0.528 ...
>    .. ..- attr(*, "dimnames")=List of 2
>    .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. ..$ : NULL
>    ..$ assign       : int [1:3] 0 1 2
>    ..$ qr           :List of 5
>    .. ..$ qr   : num [1:5, 1:3] -2.236 0.447 0.447 0.447 0.447 ...
>    .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. ..$ : chr [1:5] "1" "2" "3"
"4" ...
>    .. .. .. ..$ : chr [1:3] "(Intercept)" "lambda"
"I(lambda^2)"
>    .. .. ..- attr(*, "assign")= int [1:3] 0 1 2
>    .. ..$ qraux: num [1:3] 1.45 1.12 1.78
>    .. ..$ pivot: int [1:3] 1 2 3
>    .. ..$ tol  : num 1e-07
>    .. ..$ rank : int 3
>    .. ..- attr(*, "class")= chr "qr"
>    ..$ df.residual  : int 2
>    ..$ xlevels      : Named list()
>    ..$ call         : language lm(formula = variance.jackknife ~ lambda +
> I(lambda^2))
>    ..$ terms        :Classes 'terms', 'formula' length 3
variance.jackknife ~
> lambda + I(lambda^2)
>    .. .. ..- attr(*, "variables")= language
list(variance.jackknife, lambda,
> I(lambda^2))
>    .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
>    .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. ..$ : chr [1:3] "variance.jackknife"
"lambda" "I(lambda^2)"
>    .. .. .. .. ..$ : chr [1:2] "lambda" "I(lambda^2)"
>    .. .. ..- attr(*, "term.labels")= chr [1:2] "lambda"
"I(lambda^2)"
>    .. .. ..- attr(*, "order")= int [1:2] 1 1
>    .. .. ..- attr(*, "intercept")= int 1
>    .. .. ..- attr(*, "response")= int 1
>    .. .. ..- attr(*, ".Environment")=<environment:
0x0bd89124>
>    .. .. ..- attr(*, "predvars")= language
list(variance.jackknife, lambda,
> I(lambda^2))
>    .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.4" "numeric"
> "numeric"
>    .. .. .. ..- attr(*, "names")= chr [1:3]
"variance.jackknife" "lambda"
> "I(lambda^2)"
>    ..$ model        :'data.frame':       5 obs. of  3 variables:
>    .. ..$ variance.jackknife: num [1:5, 1:4] 1.684 1.144 0.85 0.624 0.519
...
>    .. ..$ lambda            : num [1:5] 0 0.5 1 1.5 2
>    .. ..$ I(lambda^2)       :Class 'AsIs'  num [1:5] 0 0.25 1 2.25
4
>    .. ..- attr(*, "terms")=Classes 'terms',
'formula' length 3
> variance.jackknife ~ lambda + I(lambda^2)
>    .. .. .. ..- attr(*, "variables")= language
list(variance.jackknife,
> lambda, I(lambda^2))
>    .. .. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
>    .. .. .. .. ..- attr(*, "dimnames")=List of 2
>    .. .. .. .. .. ..$ : chr [1:3] "variance.jackknife"
"lambda" "I(lambda^2)"
>    .. .. .. .. .. ..$ : chr [1:2] "lambda"
"I(lambda^2)"
>    .. .. .. ..- attr(*, "term.labels")= chr [1:2]
"lambda" "I(lambda^2)"
>    .. .. .. ..- attr(*, "order")= int [1:2] 1 1
>    .. .. .. ..- attr(*, "intercept")= int 1
>    .. .. .. ..- attr(*, "response")= int 1
>    .. .. .. ..- attr(*, ".Environment")=<environment:
0x0bd89124>
>    .. .. .. ..- attr(*, "predvars")= language
list(variance.jackknife,
> lambda, I(lambda^2))
>    .. .. .. ..- attr(*, "dataClasses")= Named chr [1:3]
"nmatrix.4" "numeric"
> "numeric"
>    .. .. .. .. ..- attr(*, "names")= chr [1:3]
"variance.jackknife" "lambda"
> "I(lambda^2)"
>    ..- attr(*, "class")= chr [1:2] "mlm" "lm"
>   $ variance.jackknife       : num [1:2, 1:2] 3.04 -0.58 -0.58 0.111
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. ..$ : chr [1:2] "(Intercept)" "x"
>   $ variance.jackknife.lambda: num [1:6, 1:5] -1 0 0.5 1 1.5 ...
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : chr [1:6] "1" "" "" ""
...
>    .. ..$ : NULL
>   $ PSI                      : num [1:1615, 1:10] -0.1032 -0.2952 -0.1847
> -0.0656 -0.3062 ...
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : chr [1:1615] "1" "2" "3"
"4" ...
>    .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
>   $ c11                      : num [1:10, 1:10] 0.151 0.78 0.144 0.748
0.139
> ...
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
>    .. ..$ : chr [1:10] "(Intercept)" "x"
"(Intercept)" "x" ...
>   $ a11                      : num [1:10, 1:10] -0.151 -0.784 0 0 0 ...
>   $ sigma                    : num [1:10, 1:10] 2979 -572 2482 -478 2132
...
>   $ sigma.gamma              : num [1:6, 1:6] 2961 -1030 192 -569 195 ...
>   $ g                        : num [1:6, 1:2] 1 -1 1 0 0 0 0 0 0 1 ...
>   $ s                        : num [1:6, 1:10] -1 0 0 0 0 0 0 0 0 -1 ...
>   $ variance.asymptotic      : num [1:2, 1:2] 3.697 -0.704 -0.704 0.134
>    ..- attr(*, "dimnames")=List of 2
>    .. ..$ : chr [1:2] "(Intercept)" "x"
>    .. ..$ : chr [1:2] "(Intercept)" "x"
>   - attr(*, "class")= chr "mcsimex"
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/How-to-Extract-Information-from-SIMEX-Output-tp3459082p3459082.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.