-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathStan-Ymet-Xnom2fac-MnormalHom-Example.R
158 lines (155 loc) · 6.88 KB
/
Stan-Ymet-Xnom2fac-MnormalHom-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
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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
# Example for Jags-Ymet-Xnom2fac-MnormalHom.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 The data file
fileNameRoot = "SalaryNormalHom-"
graphFileType = "eps"
myDataFrame = read.csv( file="Salary.csv" )
# Re-label and re-order the Pos factor:
myDataFrame$Pos = factor( myDataFrame$Pos ,
levels=c("FT3","FT2","FT1","NDW","DST") ,
ordered=TRUE ,
labels=c("Assis","Assoc","Full","Endow","Disting") )
# Specify the column names in the data file relevant to the analysis:
yName="Salary"
# x1 should be factor with fewer levels, to plot in single pane:
x1Name="Pos"
x2Name="Org"
# 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):
x1contrasts = list(
list( c("Full") , c("Assoc") , compVal=0.0 , ROPE=c(-1000,1000) ) ,
list( c("Assoc") , c("Assis") , compVal=0.0 , ROPE=c(-1000,1000) )
)
x2contrasts = list(
list( c("CHEM") , c("ENG") , compVal=0.0 , ROPE=c(-1000,1000) ) ,
list( c("CHEM") , c("PSY") , compVal=0.0 , ROPE=c(-1000,1000) ) ,
list( c("BFIN") , c("PSY","CHEM","ENG") , compVal=0.0 , ROPE=c(-1000,1000) )
)
# Each interaction contrast is a list of 2 lists of 2 vectors of level names,
# a comparison value (typically 0.0), and a ROPE (which could be NULL)::
x1x2contrasts = list(
list( list( c("Full") , c("Assis") ) ,
list( c("CHEM") , c("ENG") ) ,
compVal=0.0 , ROPE=c(-1000,1000) ) ,
list( list( c("Full") , c("Assis") ) ,
list( c("CHEM") , c("PSY") ) ,
compVal=0.0 , ROPE=c(-1000,1000) ) ,
list( list( c("Full") , c("Assoc","Assis") ) ,
list( c("BFIN") , c("PSY","CHEM","ENG") ) ,
compVal=0.0 , ROPE=c(-1000,1000) )
)
# fileNameRoot = "SplitPlotAgriData-NOSUBJ-"
# graphFileType = "eps"
# myDataFrame = read.csv( "SplitPlotAgriData.csv" )
# # Specify the column names in the data file relevant to the analysis:
# yName="Yield"
# x1Name="Fert"
# x2Name="Till"
# #xSubjectName="Field"
# x1contrasts = list(
# list( c("Deep","Surface") , c("Broad") , compVal=0.0 , ROPE=c(-5,5) )
# )
# x2contrasts = list(
# list( c("Moldbrd") , c("Ridge") , compVal=0.0 , ROPE=c(-5,5) ) ,
# list( c("Moldbrd","Ridge") , c("Chisel") , compVal=0.0 , ROPE=c(-5,5) )
# )
# x1x2contrasts = list(
# list( list( c("Broad") , c("Deep","Surface") ) ,
# list( c("Chisel","Moldbrd") , c("Ridge") ) ,
# compVal=0.0 , ROPE=c(-5,5) )
# )
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Stan-Ymet-Xnom2fac-MnormalHom.R")
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
mcmcCoda = genMCMC( datFrm=myDataFrame ,
yName=yName , x1Name=x1Name , x2Name=x2Name ,
numSavedSteps=15000 , thinSteps=5 ,
saveName=fileNameRoot )
#-------------------------------------------------------------------------------
# Display diagnostics of chain, for specified parameters:
parameterNames = varnames(mcmcCoda)
show( parameterNames ) # show all parameter names, for reference
for ( parName in c("b0","b1[1]","b2[1]","b1b2[1,1]","y_sigma") ) {
diagMCMC( codaObject=mcmcCoda , parName=parName ,
saveName=fileNameRoot , saveType=graphFileType )
}
#-------------------------------------------------------------------------------
# Get summary statistics of chain:
summaryInfo = smryMCMC( mcmcCoda ,
datFrm=myDataFrame , x1Name=x1Name , x2Name=x2Name ,
x1contrasts=x1contrasts ,
x2contrasts=x2contrasts ,
x1x2contrasts=x1x2contrasts ,
saveName=fileNameRoot )
show(summaryInfo)
# Display posterior information:
plotMCMC( mcmcCoda ,
datFrm=myDataFrame , yName=yName , x1Name=x1Name , x2Name=x2Name ,
x1contrasts=x1contrasts ,
x2contrasts=x2contrasts ,
x1x2contrasts=x1x2contrasts ,
saveName=fileNameRoot , saveType=graphFileType )
#-------------------------------------------------------------------------------
# Other specific comparisons of cells:
if ( fileNameRoot == "SalaryNormalHom-" ) {
# THIS x1level minus THAT x1level at AT x2level:
THISx1 = "Full"
THATx1 = "Assis"
ATx2 = "CHEM"
THISidx = which(levels(myDataFrame[,x1Name])==THISx1)
THATidx = which(levels(myDataFrame[,x1Name])==THATx1)
ATidx = which(levels(myDataFrame[,x2Name])==ATx2)
openGraph(height=4,width=4)
compInfo = plotPost(
as.matrix(mcmcCoda)[,paste("m[",THISidx,",",ATidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("m[",THATidx,",",ATidx,"]",sep="")] ,
main=paste(THISx1,"-",THATx1,"@",ATx2) ,
xlab=paste("Difference in",yName) ,
compVal=0 ,ROPE=c(-1000,1000) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx1,"-",THATx1,"At",ATx2,sep=""),
type=graphFileType)
# THIS x1level minus THAT x1level at AT x2level:
THISx1 = "Full"
THATx1 = "Assis"
ATx2 = "PSY"
THISidx = which(levels(myDataFrame[,x1Name])==THISx1)
THATidx = which(levels(myDataFrame[,x1Name])==THATx1)
ATidx = which(levels(myDataFrame[,x2Name])==ATx2)
openGraph(height=4,width=4)
compInfo = plotPost(
as.matrix(mcmcCoda)[,paste("m[",THISidx,",",ATidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("m[",THATidx,",",ATidx,"]",sep="")] ,
main=paste(THISx1,"-",THATx1,"@",ATx2) ,
xlab=paste("Difference in",yName) ,
compVal=0 ,ROPE=c(-1000,1000) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx1,"-",THATx1,"At",ATx2,sep=""),
type=graphFileType)
# THIS x2level minus THAT x2level at AT x1level:
THISx2 = "PSY"
THATx2 = "ENG"
ATx1 = "Full"
THISidx = which(levels(myDataFrame[,x2Name])==THISx2)
THATidx = which(levels(myDataFrame[,x2Name])==THATx2)
ATidx = which(levels(myDataFrame[,x1Name])==ATx1)
openGraph(height=4,width=4)
compInfo = plotPost(
as.matrix(mcmcCoda)[,paste("m[",ATidx,",",THISidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("m[",ATidx,",",THATidx,"]",sep="")] ,
main=paste(THISx2,"-",THATx2,"@",ATx1) ,
xlab=paste("Difference in",yName) ,
compVal=0 ,ROPE=c(-1000,1000) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx2,"-",THATx2,"At",ATx1,sep=""),
type=graphFileType)
}
#-------------------------------------------------------------------------------