DBDA2Estan/Stan-Ymet-Xnom1fac-MnormalHom-Example.R at master · jhoupt/DBDA2Estan · GitHub

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
# Example for Stan-Ymet-Xnom1fac-MnormalHom.R
#-------------------------------------------------------------------------------
# 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 The data file

# Comment out or uncomment one of the data sections below.
# In RStudio, select section then press Cntrl-Shift-C

myDataFrame = read.csv( file="FruitflyDataReduced.csv" )
# Specify the column names in the data file relevant to the analysis:
yName="Longevity"
xName="CompanionNumber"
# Specify desired contrasts.
# Each main-effect contrast is a list of 2 vectors of level names,
# a comparison value (typically 0.0), and a ROPE (which could be NULL):
contrasts = list(
  list( c("Pregnant1","Pregnant8") , c("None0") , compVal=0.0 , ROPE=c(-1.5,1.5) ) ,
  list( c("Pregnant1","Pregnant8","None0") , c("Virgin1","Virgin8") ,
        compVal=0.0 , ROPE=c(-1.5,1.5) ) ,
  list( c("Pregnant1","Pregnant8","None0") , c("Virgin1") ,
        compVal=0.0 , ROPE=c(-1.5,1.5) ) ,
  list( c("Virgin1") , c("Virgin8") , compVal=0.0 , ROPE=c(-1.5,1.5) )
)
# Specify filename root and graphical format for saving output.
# Otherwise specify as NULL or leave saveName and saveType arguments
# out of function calls.
fileNameRoot = "FruitflyData-NormalHom-"
graphFileType = "eps"

# myDataFrame = read.csv( file="NonhomogVarData.csv" )
# yName="Y"
# xName="Group"
# contrasts = list(
#   list( c("D") , c("A") , compVal=0.0 , ROPE=c(-1,1) ) ,
#   list( c("C") , c("B") , compVal=0.0 , ROPE=c(-1,1) )
# )
# fileNameRoot = "NonhomogVarData-NormalHom-"
# graphFileType = "eps"

# myDataFrame = read.csv( file="AnovaShrinkageData.csv" )
# yName="Y"
# xName="Group"
# contrasts = list(
#   list( c("U") , c("A") , compVal=0.0 , ROPE=c(-1,1) ) ,
#   list( c("M") , c("A") , compVal=0.0 , ROPE=c(-1,1) ) ,
#   list( c("G") , c("A") , compVal=0.0 , ROPE=c(-1,1) )
# )
# fileNameRoot = "AnovaShrinkageData-FixSigmaB-"
# graphFileType = "eps"

#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Stan-Ymet-Xnom1fac-MnormalHom.R")
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
mcmcCoda = genMCMC( datFrm=myDataFrame , yName=yName , xName=xName ,
                    numSavedSteps=11000 , thinSteps=1 , saveName=fileNameRoot )
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda)
show( parameterNames ) # show all parameter names, for reference
for ( parName in c("y_sigma","b0","b[1]","a_sigma") ) {
  diagMCMC( codaObject=mcmcCoda , parName=parName ,
            saveName=fileNameRoot , saveType=graphFileType )
}
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda ,
                        datFrm=myDataFrame , xName=xName ,
                        contrasts=contrasts ,
                        saveName=fileNameRoot )
show(summaryInfo)
# Display posterior information:
plotMCMC( mcmcCoda ,
          datFrm=myDataFrame , yName=yName , xName=xName ,
          contrasts=contrasts ,
          saveName=fileNameRoot , saveType=graphFileType )
#-------------------------------------------------------------------------------