sim-ar1.Rout

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
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.

> library("foreach")
> library("doParallel")
Loading required package: iterators
Loading required package: parallel
> library("parallel")
> source("init.R")

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

make: Nothing to be done for `all'.
> source("sim.R")
> 
> ## set number of Monte Carlo replicates
> M <- 1000
> 
> ## set number of threads to use for parallel processing and the random seed
> ## (nb: these two values ensure that the results are replicable)
> cores <- 2
> seed <- 0
> 
> cl <- makeCluster(getOption("cl.cores", cores))
> clusterEvalQ(cl, source("init.R"))
make: Nothing to be done for `all'.
make: Nothing to be done for `all'.
[[1]]
[[1]]$value
[[1]]$value$value
function (origin = 1, sd = 1, coralpha = sqrt(0.5), corstr = c("ar1", 
    "exchangeable"), beta0 = c(-0.2, 0, 0, 0.2, 0), beta1 = rep(0, 
    4), eta = c(0, 0, 0.8, -0.8, 0), mu = rep(0, 3), theta0 = c(0, 
    0.8), theta1 = c(0, 0), coef.avail = c(100, rep(0, 3)), coef.state = rep(0, 
    5), tfun = NULL, lag = 3 + any(beta1 != 0)) 
{
    corstr <- match.arg(corstr)
    if (is.null(tfun)) 
        tfun <- rep(list(function(tcur, tmax) rep(0, length(tcur))), 
            4)
    list(origin = 1, lag = lag, sd = sd, coralpha = coralpha, 
        corstr = corstr, beta0 = setNames(beta0, c("one", "tmod", 
            "base", "state", "lag1a")), beta1 = setNames(beta1, 
            c("one", "lag1tmod", "base", "lag1state")), eta = setNames(eta, 
            c("one", "base", "state", "lag1a", "lag1y")), mu = setNames(mu, 
            c("one", "ty", "base")), theta0 = setNames(theta0, 
            c("avail", "state")), theta1 = setNames(theta1, c("lag1avail", 
            "lag1state")), coef.avail = setNames(coef.avail, 
            c("one", "tavail", "lag1a", "lag1y")), coef.state = setNames(coef.state, 
            c("one", "tstate", "base", "lag1state", "lag1a")), 
        tfun = setNames(tfun, c("ty", "tmod", "tavail", "tstate")))
}

[[1]]$value$visible
[1] FALSE


[[1]]$visible
[1] FALSE


[[2]]
[[2]]$value
[[2]]$value$value
function (origin = 1, sd = 1, coralpha = sqrt(0.5), corstr = c("ar1", 
    "exchangeable"), beta0 = c(-0.2, 0, 0, 0.2, 0), beta1 = rep(0, 
    4), eta = c(0, 0, 0.8, -0.8, 0), mu = rep(0, 3), theta0 = c(0, 
    0.8), theta1 = c(0, 0), coef.avail = c(100, rep(0, 3)), coef.state = rep(0, 
    5), tfun = NULL, lag = 3 + any(beta1 != 0)) 
{
    corstr <- match.arg(corstr)
    if (is.null(tfun)) 
        tfun <- rep(list(function(tcur, tmax) rep(0, length(tcur))), 
            4)
    list(origin = 1, lag = lag, sd = sd, coralpha = coralpha, 
        corstr = corstr, beta0 = setNames(beta0, c("one", "tmod", 
            "base", "state", "lag1a")), beta1 = setNames(beta1, 
            c("one", "lag1tmod", "base", "lag1state")), eta = setNames(eta, 
            c("one", "base", "state", "lag1a", "lag1y")), mu = setNames(mu, 
            c("one", "ty", "base")), theta0 = setNames(theta0, 
            c("avail", "state")), theta1 = setNames(theta1, c("lag1avail", 
            "lag1state")), coef.avail = setNames(coef.avail, 
            c("one", "tavail", "lag1a", "lag1y")), coef.state = setNames(coef.state, 
            c("one", "tstate", "base", "lag1state", "lag1a")), 
        tfun = setNames(tfun, c("ty", "tmod", "tavail", "tstate")))
}

[[2]]$value$visible
[1] FALSE


[[2]]$visible
[1] FALSE


> registerDoParallel(cl)
> 
> ## control parameters for the generative model
> beta0 <- c(-0.2, 0, 0, 0, 0)     # unmoderated proximal effect
> beta1 <- c(-0.1, 0, 0, 0)        # unmoderated delayed effect
> coef.state <- c(0, 0, 0, 0, 0.1) # state depends on past treatment
> eta <- rep(0, 5)                 # treatment probability = 1/2
> 
> sim.ar1 <- function() {
+   out <- NULL
+   for (n in c(30, 60)) {
+     for (tmax in c(30, 50)) {
+       ## obtain true correlation matrix, trimmed down to effective size
+       ## ("effective" observations avoid (lags of) initial values)
+       attrib <- attributes(rsnmm(n, tmax, beta0 = beta0, beta1 = beta1,
+                                  coef.state = coef.state, eta = eta))
+       ## by default the true correlation structure is AR(1) with (u, t)th error
+       ## correlation sqrt(0.5)^abs(u - t)
+       cormatrix <- attrib$cormatrix[1:(tmax - attrib$lag + 1),
+                                     1:(tmax - attrib$lag + 1)]
+       clusterSetRNGStream(cl, seed)
+       out <- rbind(out,
+                    sim(n, tmax, M,
+                        ## regress response on proximal treatment, centered by the
+                        ## true treatment probability
+                        y.formula = list(indep = y ~ state + I(a - prob),
+                                         ar1 = y ~ state + I(a - prob)),
+                        y.names = c(indep = "Independence",
+                                    ar1 = "AR(1)"),
+                        y.label = list(indep = "I(a - prob)",
+                                       ar1 = "I(a - prob)"),
+                        ## employ different working correlation structures
+                        y.args = list(indep = list(),
+                                      ar1 = list(corstr = "userdefined",
+                                                 wcor = cormatrix)),
+                        a.formula = NULL, a.names = NULL,
+                        beta0 = beta0, beta1 = beta1, coef.state = coef.state,
+                        eta = eta))
+     }
+   }
+   out
+ }
> 
> ar1 <- sim.ar1()

Generative model attributes

$origin
[1] 1

$lag
[1] 4

$sd
[1] 1

$coralpha
[1] 0.7071068

$corstr
[1] "ar1"

$beta0
  one  tmod  base state lag1a 
 -0.2   0.0   0.0   0.0   0.0 

$beta1
      one  lag1tmod      base lag1state 
     -0.1       0.0       0.0       0.0 

$eta
  one  base state lag1a lag1y 
    0     0     0     0     0 

$mu
 one   ty base 
   0    0    0 

$theta0
avail state 
  0.0   0.8 

$theta1
lag1avail lag1state 
        0         0 

$coef.avail
   one tavail  lag1a  lag1y 
   100      0      0      0 

$coef.state
      one    tstate      base lag1state     lag1a 
      0.0       0.0       0.0       0.0       0.1 

$tfun
$tfun$ty
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x33d66a8>

$tfun$tmod
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x33d66a8>

$tfun$tavail
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x33d66a8>

$tfun$tstate
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x33d66a8>


Analysis models

  Independence: y ~ state + I(a - prob)

  AR(1): y ~ state + I(a - prob)

Treatment probability models


Generative model attributes

$origin
[1] 1

$lag
[1] 4

$sd
[1] 1

$coralpha
[1] 0.7071068

$corstr
[1] "ar1"

$beta0
  one  tmod  base state lag1a 
 -0.2   0.0   0.0   0.0   0.0 

$beta1
      one  lag1tmod      base lag1state 
     -0.1       0.0       0.0       0.0 

$eta
  one  base state lag1a lag1y 
    0     0     0     0     0 

$mu
 one   ty base 
   0    0    0 

$theta0
avail state 
  0.0   0.8 

$theta1
lag1avail lag1state 
        0         0 

$coef.avail
   one tavail  lag1a  lag1y 
   100      0      0      0 

$coef.state
      one    tstate      base lag1state     lag1a 
      0.0       0.0       0.0       0.0       0.1 

$tfun
$tfun$ty
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x2db23d0>

$tfun$tmod
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x2db23d0>

$tfun$tavail
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x2db23d0>

$tfun$tstate
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x2db23d0>


Analysis models

  Independence: y ~ state + I(a - prob)

  AR(1): y ~ state + I(a - prob)

Treatment probability models


Generative model attributes

$origin
[1] 1

$lag
[1] 4

$sd
[1] 1

$coralpha
[1] 0.7071068

$corstr
[1] "ar1"

$beta0
  one  tmod  base state lag1a 
 -0.2   0.0   0.0   0.0   0.0 

$beta1
      one  lag1tmod      base lag1state 
     -0.1       0.0       0.0       0.0 

$eta
  one  base state lag1a lag1y 
    0     0     0     0     0 

$mu
 one   ty base 
   0    0    0 

$theta0
avail state 
  0.0   0.8 

$theta1
lag1avail lag1state 
        0         0 

$coef.avail
   one tavail  lag1a  lag1y 
   100      0      0      0 

$coef.state
      one    tstate      base lag1state     lag1a 
      0.0       0.0       0.0       0.0       0.1 

$tfun
$tfun$ty
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3342868>

$tfun$tmod
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3342868>

$tfun$tavail
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3342868>

$tfun$tstate
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3342868>


Analysis models

  Independence: y ~ state + I(a - prob)

  AR(1): y ~ state + I(a - prob)

Treatment probability models


Generative model attributes

$origin
[1] 1

$lag
[1] 4

$sd
[1] 1

$coralpha
[1] 0.7071068

$corstr
[1] "ar1"

$beta0
  one  tmod  base state lag1a 
 -0.2   0.0   0.0   0.0   0.0 

$beta1
      one  lag1tmod      base lag1state 
     -0.1       0.0       0.0       0.0 

$eta
  one  base state lag1a lag1y 
    0     0     0     0     0 

$mu
 one   ty base 
   0    0    0 

$theta0
avail state 
  0.0   0.8 

$theta1
lag1avail lag1state 
        0         0 

$coef.avail
   one tavail  lag1a  lag1y 
   100      0      0      0 

$coef.state
      one    tstate      base lag1state     lag1a 
      0.0       0.0       0.0       0.0       0.1 

$tfun
$tfun$ty
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3372838>

$tfun$tmod
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3372838>

$tfun$tavail
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3372838>

$tfun$tstate
function (tcur, tmax) 
rep(0, length(tcur))
<environment: 0x3372838>


Analysis models

  Independence: y ~ state + I(a - prob)

  AR(1): y ~ state + I(a - prob)

Treatment probability models

> save(ar1, file = "sim-ar1.RData")
> 
> stopCluster(cl)
> 
> proc.time()
     user    system   elapsed 
    7.160     0.287 12534.663