-
Notifications
You must be signed in to change notification settings - Fork 14
/
Stan-Ycount-Xnom2fac-MpoissonExp-Example.R
187 lines (187 loc) · 7.86 KB
/
Stan-Ycount-Xnom2fac-MpoissonExp-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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
# 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
#.............................................................................
myDataFrame = read.csv( file="HairEyeColor.csv" )
# Alter count by a multiplier, for demo purposes:
countMult = 1
myDataFrame$Count = round( myDataFrame$Count * countMult )
fileNameRoot = paste0("HairEyeColor-",countMult,"-")
yName="Count"
x1Name="Eye"
x2Name="Hair"
x1contrasts = list(
list( c("Green") , c("Hazel") , compVal=0.0 , ROPE=c(-0.1,0.1) ) ,
list( c("Blue") , c("Green") , compVal=0.0 , ROPE=c(-0.1,0.1) )
)
x2contrasts = list(
list( c("Black") , c("Blond") , compVal=0.0 , ROPE=c(-0.1,0.1) ) ,
list( c("Brown") , c("Red") , compVal=0.0 , ROPE=c(-0.1,0.1) )
)
x1x2contrasts = list(
list( list( c("Blue") , c("Brown") ) ,
list( c("Black") , c("Blond") ) ,
compVal=0.0 , ROPE=c(-0.1,0.1) )
)
numSavedSteps = 1200
thinSteps = 1
#.............................................................................
# myDataFrame = read.csv( file="FourByFourCount.csv" )
# # Increase count by an integer multiplier:
# countMult = 11
# myDataFrame$Count = myDataFrame$Count * countMult
# fileNameRoot = paste0("FourByFourCount-",countMult,"-")
# yName="Count"
# x1Name="A"
# x2Name="B"
# x1contrasts = list(
# list( c("A1","A2") , c("A3","A4") , compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# x2contrasts = list(
# list( c("B1","B2") , c("B3","B4") , compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# x1x2contrasts = list(
# list( list( c("A1","A2") , c("A3","A4") ) ,
# list( c("B1","B2") , c("B3","B4") ) ,
# compVal=0.0 , ROPE=c(-0.1,0.1) ) ,
# list( list( c("A2") , c("A3") ) ,
# list( c("B2") , c("B3") ) ,
# compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# numSavedSteps = 12000
# thinSteps = 10
#.............................................................................
# myDataFrame = read.csv( file="CrimeDrink.csv" )
# # Alter count by a multiplier, for demo purposes:
# countMult = 1
# myDataFrame$Count = round( myDataFrame$Count * countMult )
# fileNameRoot = paste0("CrimeDrink-",countMult,"-")
# yName="Count"
# x2Name="Drink"
# x1Name="Crime"
# x2contrasts = list(
# list( c("Drinker") , c("Nondrink") , compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# x1contrasts = list(
# list( c("Fraud") , c("Violence") , compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# x1x2contrasts = list(
# list( list( c("Fraud") , c("Violence") ) ,
# list( c("Drinker") , c("Nondrink") ) ,
# compVal=0.0 , ROPE=c(-0.1,0.1) )
# )
# numSavedSteps = 12000
# thinSteps = 10
#.............................................................................
graphFileType = "eps"
#-------------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Stan-Ycount-Xnom2fac-MpoissonExp.R")
#-------------------------------------------------------------------------------
# Generate the MCMC chain:
mcmcCoda = genMCMC( datFrm=myDataFrame ,
yName=yName , x1Name=x1Name , x2Name=x2Name ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
saveName=fileNameRoot, nChains=4 )
#-------------------------------------------------------------------------------
# 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]","pp_x1x2_p[1,1]",
"a1_sd","a1a2_sd") ) {
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 ( substr(fileNameRoot,1,12) == "HairEyeColor" ) {
# THIS x1level minus THAT x1level at AT x2level:
THISx1 = "Blue"
THATx1 = "Brown"
ATx2 = "Black"
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("b1b2[",THISidx,",",ATidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("b1b2[",THATidx,",",ATidx,"]",sep="")] ,
main=paste(THISx1,"-",THATx1,"@",ATx2) ,
xlab=paste("Beta Deflect. Diff.") ,
compVal=0 ,ROPE=c(-0.1,0.1) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx1,"-",THATx1,"At",ATx2,sep=""),
type=graphFileType)
# THIS x1level minus THAT x1level at AT x2level:
THISx1 = "Blue"
THATx1 = "Brown"
ATx2 = "Blond"
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("b1b2[",THISidx,",",ATidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("b1b2[",THATidx,",",ATidx,"]",sep="")] ,
main=paste(THISx1,"-",THATx1,"@",ATx2) ,
xlab=paste("Beta Deflect. Diff.") ,
compVal=0 ,ROPE=c(-0.1,0.1) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx1,"-",THATx1,"At",ATx2,sep=""),
type=graphFileType)
# THIS x2level minus THAT x2level at AT x1level:
THISx2 = "Black"
THATx2 = "Blond"
ATx1 = "Blue"
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("b1b2[",ATidx,",",THISidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("b1b2[",ATidx,",",THATidx,"]",sep="")] ,
main=paste(THISx2,"-",THATx2,"@",ATx1) ,
xlab=paste("Beta Deflect. Diff.") ,
compVal=0 ,ROPE=c(-0.1,0.1) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx2,"-",THATx2,"At",ATx1,sep=""),
type=graphFileType)
# THIS x2level minus THAT x2level at AT x1level:
THISx2 = "Black"
THATx2 = "Blond"
ATx1 = "Brown"
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("b1b2[",ATidx,",",THISidx,"]",sep="")] -
as.matrix(mcmcCoda)[,paste("b1b2[",ATidx,",",THATidx,"]",sep="")] ,
main=paste(THISx2,"-",THATx2,"@",ATx1) ,
xlab=paste("Beta Deflect. Diff.") ,
compVal=0 ,ROPE=c(-0.1,0.1) )
show(compInfo)
saveGraph(file=paste(fileNameRoot,THISx2,"-",THATx2,"At",ATx1,sep=""),
type=graphFileType)
}
#-------------------------------------------------------------------------------