-
Notifications
You must be signed in to change notification settings - Fork 14
/
Stan-Ymet-XmetMulti-MrobustShrink-Example.R
executable file
·81 lines (79 loc) · 3.85 KB
/
Stan-Ymet-XmetMulti-MrobustShrink-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
# Example for Stan-Ymet-XmetMulti-MrobustShrink.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 data file; specify column names of x (predictors) and y (predicted):
myData = read.csv( file="Guber1999data.csv" )
# UNCOMMENT ONE OF THE FOLLOWING SECTIONS (In RStudio, select and ctrl-shift-C)
# #.........................................................................
# # Two predictors:
# yName = "SATT" ; xName = c("Spend","PrcntTake")
# fileNameRoot = "Guber1999data-Stan-Shrink-"
# numSavedSteps=15000 ; thinSteps=5
#...........................................................................
# # Two predictors with redundant predictor:
# PropNotTake = (100-myData[,"PrcntTake"])/100
# myData = cbind( myData , PropNotTake )
# yName = "SATT" ; xName = c("Spend","PrcntTake","PropNotTake")
# fileNameRoot = "Guber1999data-Stan-Shrink-Redund-"
# numSavedSteps=15000 ; thinSteps=15
#...........................................................................
# # Two predictors with two redundant predictors:
# PropNotTake = (100-myData[,"PrcntTake"])/100
# Partic = myData[,"PrcntTake"]/10
# myData = cbind( myData , PropNotTake , Partic )
# yName = "SATT" ; xName = c("Spend","PrcntTake","PropNotTake","Partic")
# fileNameRoot = "Guber1999data-Stan-Shrink-Redund2-"
# numSavedSteps=15000 ; thinSteps=15
#...........................................................................
# # Four predictors:
# yName = "SATT" ; xName = c("Spend","PrcntTake","StuTeaRat","Salary")
# fileNameRoot = "Guber1999data-Stan-Shrink-4X-"
# numSavedSteps=15000 ; thinSteps=20
#.............................................................................
# Append columns of random predictors:
standardize = function(v){(v-mean(v))/sd(v)}
Ny=nrow(myData)
NxRand = 12
set.seed(47405)
for ( xIdx in 1:NxRand ) {
xRand = standardize(rnorm(Ny))
myData = cbind( myData , xRand )
colnames(myData)[ncol(myData)] = paste0("xRand",xIdx)
}
yName = "SATT" ; xName = c("Spend","PrcntTake", paste0("xRand",1:NxRand) )
fileNameRoot = "Guber1999data-Stan-Shrink-RandX-"
numSavedSteps=15000 ; thinSteps=5
#...........................................................................
graphFileType = "eps"
#---------------------------------------------------------------------------
# Load the relevant model into R's working memory:
source("Stan-Ymet-XmetMulti-MrobustShrink.R")
#---------------------------------------------------------------------------
# Generate the MCMC chain:
#startTime = proc.time()
mcmcCoda = genMCMC( data=myData , xName=xName , yName=yName ,
numSavedSteps=numSavedSteps , thinSteps=thinSteps ,
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 parameterNames ) {
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 ,
pairsPlot=TRUE , showCurve=FALSE ,
saveName=fileNameRoot , saveType=graphFileType )
#---------------------------------------------------------------------------