## Simulation of the variable slopes model.

## Structural Parameters

I <- 100                                #Number of students
Tmax <- 10                              #Maximum time periods

## Experimental Design:  Time and dosage patterns.
## For the moment assume all have the same number of measurements.
Time <- matrix(1/Tmax,I,Tmax-1)         #Equally spaced measurements
Dose <- matrix(0,I,Tmax-1)              #Start with 0, we will fill
                                        #this in later.
T <- rep(Tmax,I)                        #Number of measurement points
                                        #for each person
## Place holder for data
Obs <- matrix(NA,I,Tmax)

## Hyperparameters
T0 <- 1+1/Tmax                          #Effective start time
slope_mu_mu <- 1                        #Mean slope
slope_mu_std <- .1
T0k <- slope_mu_mu*T0                      #Offset to make starting mean 0

## Evidence Model priors
obs_rel <- .8
obs_mean_t1 <- 50
obs_std_t1 <- 10
res_std <- obs_std_t1*sqrt(1-obs_rel)
obs_slope <- obs_std_t1*sqrt(obs_rel)
obs_int <- obs_mean_t1

## Parameters:  Brownian Motion
var_innov <- .1                          #Variance of the innovations
treat_eff <- .25                        #Treatment effect
T0_std_err <- .25                       #Std error at time 0
slope_std <- .33                        #Std of slope

## Level 2 parameters
slope_mu <- rnorm(1,slope_mu_mu,slope_mu_std)

## Level 1 parameters
slope <- rnorm(I,slope_mu,slope_std)
T0_err <- rnorm(I,0,T0_std_err)
innov <- matrix(NA,I,Tmax-1)

## Set up dosage based on slope.
pdose <- ifelse(slope>=slope_mu,0,1-2*pnorm(slope,slope_mu,slope_std))
for (i in 1:I) {
  Dose[i,] <- rbinom(T[i]-1,1,pdose[i])*Time[i,]
}

theta <- matrix(NA,I,Tmax)
theta[,1] <- slope*T0 - T0k + T0_err
for (t in 2:Tmax) {
  innov[,t-1] <- rnorm(I,0,sqrt(var_innov*Time[,t-1]))
  theta[,t] <- theta[,t-1] + slope*Time[,t-1] +
    treat_eff*Dose[,t-1]+innov[,t-1]

}

for (t in 1:Tmax) {
  Obs[,t] <- rnorm(I,obs_int+obs.slope*theta[,t],res.std)
}