-
Notifications
You must be signed in to change notification settings - Fork 14
/
Stan-Ymet-XmetSsubj-MrobustHier-Example.R
52 lines (48 loc) · 2.51 KB
/
Stan-Ymet-XmetSsubj-MrobustHier-Example.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# Example for Stan-Ymet-XmetSsubj-MrobustHier.R
# Adapted for Stan by Joe Houpt
#-------------------------------------------------------------------------------
# Optional generic preliminaries:
graphics.off() # This closes all of R's graphics windows.
rm(list=ls()) # Careful! This clears all of R's memory!
#-------------------------------------------------------------------------------
# Load data file and specity column names of x (predictor) and y (predicted):
myData = read.csv( file="IncomeFamszState3yr.csv" , comment.char="#")
xName = "FamilySize" ; yName = "MedianIncome" ; sName="State" ; wName="SampErr"
fileNameRoot = "IncomeFamszState3yr-Stan-"
# myData = read.csv( file="HierLinRegressData.csv" )
# xName = "X" ; yName = "Y" ; sName="Subj" ; wName=NULL
# fileNameRoot = "HierLinRegressData-Stan-"
# myData = read.csv( file="BugsRatsData.csv" )
# xName = "Day" ; yName = "Weight" ; sName="Subj" ; wName=NULL
# fileNameRoot = "BugsRatsData-Stan-"
graphFileType = "eps"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Stan-Ymet-XmetSsubj-MrobustHier.R")
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
#startTime = proc.time()
mcmcCoda = genMCMC( data=myData ,
xName=xName , yName=yName , sName=sName , wName=wName ,
numSavedSteps=12000 , thinSteps=5 , saveName=fileNameRoot )
#stopTime = proc.time()
#duration = stopTime - startTime
#show(duration)
# #-------------------------------------------------------------------------------
# # Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda) # get all parameter names
for ( parName in c("beta_0_mu","beta_1_mu","nu","sigma",
"beta_0[1]","beta_1[1]") ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda , saveName=fileNameRoot )
show(summaryInfo)
# Display posterior information:
plotMCMC( mcmcCoda , data=myData ,
xName=xName , yName=yName , sName=sName , wName=wName ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
#-------------------------------------------------------------------------------