## Simulation of the variable slopes model, independent errors

## 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,Tmax-1,I)         #Equally spaced measurements
## vector[I] Time[Tmax-1];
Dose <- matrix(0,Tmax-1,I)              #Start with 0, we will fill
## vector[I] Dose[Tmax-1];
                                        #this in later.
T <- rep(Tmax,I)                        #Number of measurement points
                                        #for each person
## Place holder for data
ObsA <- matrix(NA,Tmax,I)
## vector[I] Obs[Tmax-1];

## Hyperparameters
slope_mu_mu <- 1                        #Mean slope
slope_mu_std <- .1

## 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
slope_std <- .33
slope_r2 <- .36
## Level 2 parameters
slope_mu <- rnorm(1,slope_mu_mu,slope_mu_std)

## Level 1 parameters
slope_res <- rnorm(I,0.0,1.0)
T0_val <- rnorm(I,0,1.0)
innovA <- matrix(NA,Tmax-1,I)
## vector[I] innov[Tmax-1];
slope <- slope_mu +
  slope_std*(sqrt(1-slope_r2)*slope_res + slope_r2*T0_val)


## 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]
}

cumTime <- rbind(0,apply(Time,2,cumsum))
cumDose <- rbind(0,apply(Dose,2,cumsum))


thetaA <- matrix(NA,Tmax,I)
## vector[I] theta[Tmax];
thetaA[1,] <- T0_val
for (t in 2:Tmax) {
  innovA[t-1,] <- rnorm(I,0,sqrt(var_innov))
  thetaA[t,] <- T0_val + slope*cumTime[t-1,] +
    treat_eff*cumDose[t-1,]+innovA[t-1,]

}

for (t in 1:Tmax) {
  ObsA[t,] <- rnorm(I,obs_int+obs_slope*thetaA[t,],res_std)
}