########## R script: ratsModel3.Rs ########## # For fitting a random intercepts and slopes linear mixed model # to the `rats' data, with a quadratic fixed effects # component. # Last changed: 01 AUG 2024 # Set flags: doMCMC <- TRUE ; compileCode <- TRUE doSummMCMC <- TRUE doFitPlot <- TRUE doResidPlot <- TRUE # 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: nBurn <- 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] u[m]; cov_matrix[2] Sigma; real sigmaEps; vector[2] aSigma; } transformed parameters { vector[numObs] meanVec; cov_matrix[2] SigmaInv; meanVec = X*beta + to_vector(u[idnum,1]) + to_vector(u[idnum,2]).*col(X,2); SigmaInv = inverse(Sigma); } model { vector[2] scaleSigma; y ~ normal(meanVec,sigmaEps); u ~ 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); }' # 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=nBurn, iter=(nBurn + nKept),chains=1,thin=nThin,refresh=200, fit=stanCompilObj) # Save and extract relevant MCMC samples: beta0MCMC <- as.vector(extract(stanObj,"beta[1]",permuted=FALSE)) beta1MCMC <- as.vector(extract(stanObj,"beta[2]",permuted=FALSE)) beta2MCMC <- as.vector(extract(stanObj,"beta[3]",permuted=FALSE)) sigmaU0MCMC <- sqrt(as.vector(extract(stanObj,"Sigma[1,1]",permuted=FALSE))) sigmaU1MCMC <- sqrt(as.vector(extract(stanObj,"Sigma[2,2]",permuted=FALSE))) rhoMCMC <- as.vector(extract(stanObj,"Sigma[1,2]",permuted=FALSE))/(sigmaU0MCMC*sigmaU1MCMC) sigmaEpsMCMC <- as.vector(extract(stanObj,"sigmaEps",permuted=FALSE)) u0MCMC <- NULL ; u1MCMC <- NULL for (i in 1:m) { charVar <- paste("u[",as.character(i),",1]",sep = "") u0MCMC <- rbind(u0MCMC,extract(stanObj,charVar,permuted = FALSE)) charVar <- paste("u[",as.character(i),",2]",sep = "") u1MCMC <- rbind(u1MCMC,extract(stanObj,charVar,permuted = FALSE)) } # Convert to original units: beta1OrigMCMC <- beta1MCMC*(sd.y/sd.x) beta2OrigMCMC <- beta2MCMC*(sd.y/sd.x^2) beta0OrigMCMC <- mean.y + sd.y*beta0MCMC sigmaU0OrigMCMC <- sigmaU0MCMC*sd.y sigmaU1OrigMCMC <- sigmaU1MCMC*sd.y sigmaEpsOrigMCMC <- sigmaEpsMCMC*sd.y } if (doSummMCMC) { parNamesVal <- list(c(expression(beta[0])),c(expression(beta[1])),c(expression(beta[2])), c(expression(sigma[u0])),c(expression(sigma[u1])),c(expression(rho)), c(expression(sigma[epsilon]))) MCMClist <- list(cbind(beta0OrigMCMC,beta1OrigMCMC,beta2OrigMCMC, sigmaU0OrigMCMC,sigmaU1OrigMCMC,rhoMCMC, sigmaEpsOrigMCMC)) summMCMC(MCMClist,parNames=parNamesVal) readline("Hit Enter to continue.\n") } if (doFitPlot) { library(lattice) ng <- 101 xg <- seq(min(x),max(x),length=ng) fitg <- vector("list",m) credLowg <- vector("list",m) credUppg <- vector("list",m) for (i in 1:m) { fitMCMC <- (outer(rep(1,ng),(beta0MCMC + u0MCMC[i,])) + outer(xg,(beta1MCMC + u1MCMC[i,])) + outer((xg^2),beta2MCMC)) fitg[[i]] <- apply(fitMCMC,1,mean) credLowg[[i]] <- apply(fitMCMC,1,quantile,0.025) credUppg[[i]] <- apply(fitMCMC,1,quantile,0.975) } # Convert to original units (note that mean(x) is zero for these # data, which makes the next commands simpler than for the # general situation): xgOrig <- sd.x*xg fitgOrig <- vector("list",m) credLowgOrig <- vector("list",m) credUppgOrig <- vector("list",m) for (i in 1:m) { fitgOrig[[i]] <- mean.y + sd.y*fitg[[i]] credLowgOrig[[i]] <- mean.y + sd.y*credLowg[[i]] credUppgOrig[[i]] <- mean.y + sd.y*credUppg[[i]] } fitPlot <- xyplot(yOrig~xOrig|as.factor(idnum),groups=idnum, data=data.frame(x=x,y=y,idnum=idnum), layout=layoutVec,main="",strip=FALSE, xlab=xlabVal,ylab=ylabVal,as.table=TRUE, panel=function(x,y,subscripts,groups) { i <- panel.number() panel.grid() panel.polygon(c(xgOrig,rev(xgOrig)),c(credLowgOrig[[i]],rev(credUppgOrig[[i]])), col="palegreen",border=FALSE) panel.superpose(x[order(x)],y[order(x)],subscripts,groups, col="blue",cex=0.5) panel.xyplot(xgOrig,fitgOrig[[i]],col="darkgreen",type="l",lwd=2) }) print(fitPlot) readline("Hit Enter to continue.\n") } if (doResidPlot) { fittedVec <- NULL for (i in 1:m) { # Obtain current x and y vectors: xCurr <- x[idnum==i] yCurr <- y[idnum==i] # Obtain current fitted value vector: fitCurrMCMC <- (outer(rep(1,length(xCurr)),(beta0MCMC + u0MCMC[i,])) + outer(xCurr,(beta1MCMC + u1MCMC[i,])) + outer((xCurr^2),beta2MCMC)) fitCurr <- apply(fitCurrMCMC,1,mean) fittedVec <- c(fittedVec,fitCurr) } residVec <- y - fittedVec fittedVecOrig <- mean.y + sd.y*fittedVec residVecOrig <- sd.y*residVec plot(fittedVecOrig,residVecOrig,bty="l",xlab="fitted values",ylab="residuals", col="blue",) abline(0,0,col="red",lwd=2) } ######### End of ratsModel3.Rs ##########