########## R script: ratsModel3LOOana ########## # For fitting Model 3 to the rats data and obtaining # the leave-one-out object. # Last changed: 20 SEP 2022 # Set flags: doMCMC <- TRUE ; compileCode <- TRUE doLOOcalc <- TRUE options(warn=-1) # Assign data: m <- 30 uniqWeeks <- c(-14,-7,0,7,14) weeks <- rep(uniqWeeks,m) idnum <- rep(1:m,each=5) weight <- c(151,199,246,283,320, 145,199,249,293,354, 147,214,263,312,328, 155,200,237,272,297, 135,188,230,280,323, 159,210,252,298,331, 141,189,231,275,305, 159,201,248,297,338, 177,236,285,350,376, 134,182,220,260,296, 160,208,261,313,352, 143,188,220,273,314, 154,200,244,289,325, 171,221,270,326,358, 163,216,242,281,312, 160,207,248,288,324, 142,187,234,280,316, 156,203,243,283,317, 157,212,259,307,336, 152,203,246,286,321, 154,205,253,298,334, 139,190,225,267,302, 146,191,229,272,302, 157,211,250,285,323, 132,185,237,286,331, 160,207,257,303,345, 169,216,261,295,333, 157,205,248,289,316, 137,180,219,258,291, 153,200,244,286,324) xOrig <- weeks yOrig <- weight xlabVal <- "number of weeks" ylabVal <- "weight" layoutVec <- c(6,5) # Standardise for Bayesian analysis (note that # the predictor data has a sample mean of zero): sd.x <- sd(xOrig) mean.y <- mean(yOrig) ; sd.y <- sd(yOrig) x <- xOrig/sd.x ; y <- (yOrig - mean.y)/sd.y # Obtain X matrix: X <- cbind(1,x,x^2) # Set dimension values: numObs <- length(y) m <- length(unique(idnum)) # Do Markov chain Monte Carlo if requested: if (doMCMC) { library(rstan) # Set MCMC sample size paramaters: nWarm <- 1000 nKept <- 1000 nThin <- 1 ratsModel3 <- 'data { int numObs; int m; matrix[numObs,3] X; real y[numObs]; int idnum[numObs]; } transformed data { vector[2] zeroVec; zeroVec = rep_vector(0,2); } parameters { vector[3] beta; vector[2] uLin[m]; cov_matrix[2] Sigma; real sigmaEps; vector[2] aSigma; } transformed parameters { vector[numObs] meanVec; cov_matrix[2] SigmaInv; meanVec = X*beta + to_vector(uLin[idnum,1]) + to_vector(uLin[idnum,2]).*col(X,2); SigmaInv = inverse(Sigma); } model { vector[2] scaleSigma; y ~ normal(meanVec,sigmaEps); uLin ~ multi_normal(zeroVec,Sigma); beta ~ normal(0,100000); aSigma ~ inv_gamma(0.5,0.0000000001); for (k in 1:2) scaleSigma[k] = 4/aSigma[k]; Sigma ~ inv_wishart(3,diag_matrix(scaleSigma)); sigmaEps ~ cauchy(0,10000); } generated quantities { vector[numObs] log_lik; for (iObs in 1:numObs) { log_lik[iObs] = normal_lpdf(y[iObs]|X[iObs]*beta + uLin[idnum[iObs],1] + X[iObs,2]*uLin[idnum[iObs],2],sigmaEps); } } ' # Set up data list: allData <- list(y=y,X=X,idnum=idnum,numObs=numObs,m=m) # Compile code for model if required: if (compileCode) stanCompilObj <- stan(model_code=ratsModel3,data=allData, iter=1,chains=1) # Perform MCMC: stanObj <- stan(model_code=ratsModel3,data=allData,warmup=nWarm, iter=(nWarm + nKept),chains=1,thin=nThin,refresh=200, fit=stanCompilObj) } if (doLOOcalc) { library(loo) LOOobj <- loo(extract_log_lik(stanObj)) cat("\n\n\n The leave-one-out information criterion for Model 3 is:",LOOobj$estimates[3,1]) cat("\n\n\n") } ######### End of ratsModel3LOOana ##########