MAJUG Fall Meeting, Thursday, Sept. 28, 2017, SAS Office, 1530 Wilson Blvd., from 8:30am-12:30pm

[#1]

Mid-Atlantic JMP® Users Group (MAJUG) Fall Meeting Location: SAS Office, 1530 Wilson Blvd. Suite 800, Arlington-Rosslyn, VA Thursday, September 28, 2017 from 8:30-12:30 AGENDA*:

8:30 Registration/Sign-in 8:45-9:00 Welcome and Introduction  9:00-10:00 Data Transformation (Tom Donnelly, JMP Systems Engineer)  10:00-10:15 Break 10:15-10:30 MAJUG Feedback – survey, future meetings, topics, user involvement, JMP community, meeting location

10:30-11:15 Modeling Consumer Expenditure Survey Interview Data (Josh Klick, US Bureau of Labor Statistics and MAJUG Co-chair)

11:15-11:30 Break

11:30-12:00 JMP® and R Integration Application using Alaskan Pipeline data and R code from http://www.itl.nist.gov/div898/handbook/pmd/section6/pmd621.htm  (Mel Alexander, Social Security Administration and MAJUG Co-chair)

12:00-12:25 Q&A, Discussion, and Feedback 12:30 Adjourn * This schedule is subject to change without notice

MAJUG meeting will be located at the SAS Office, 1530 Wilson Blvd. (Suite 800), Arlington-Rosslyn, VA

We encourage you to attend in person to get full advantage of the meeting and learning from other JMP users. If you are unable to attend you can Register to attend using WebEx . 

Please reply to Joshua Klick at klick.joshua@bls.gov   

by 12pm Monday September 25th indicating whether you plan to attend in person or using WebEx. Instructions to accessWebEx session will be sent to WebEx registrants. 

Bios

Tom Donnelly works as a Systems Engineer in the JMP® division of SAS®. A physicist by training, he has been actively using and teaching Design of Experiments (DOE) methods for the past 25 years to develop and optimize products, processes and technologies.

Donnelly joined JMP after working as an analyst for the Modeling, Simulation & Analysis Branch of the US Army’s Edgewood Chemical Biological Center at Aberdeen Proving Ground, MD. There, he used DOE to develop, test, and evaluate technologies for detection, protection and decontamination of chemical and biological agents.

Josh Klick is a Senior Economist with the Bureau of Labor Statistics in the Division of Consumer Prices and Price Indices. His work focuses on the methodology and systems development that supports the construction of the Consumer Price Index (CPI). He also co-chairs the Mid-Atlantic JMP® Users Group (MAJUG). He earned a master’s degree in Public Policy- Economic Development from the University of Minnesota-Twin Cities. He has certifications as a SAS Basic and Advanced Programmer.

Mel Alexander is an Operations Research Analyst with the Social Security Administration (SSA) and Analytician at the University of Maryland’s Medical Center (UMMC). At SSA, he uses SAS® software to ensure proper payments go to eligible beneficiaries. At UMMC, he uses JMP® statistical analytics to help radiologists and neurosurgeons identify key factors that increase patient’s chances of surviving traumatic head, spinal cord, and abdominal injuries. He co-chairs MAJUG and is a section co-chair for the South East SAS Users Group (SESUG) 2017 Conference. He earned a master’s degree in Biostatistics from the University of North Carolina. He is an American Statistical Association member, an American Society for Quality (ASQ) Fellow and Certified Quality Engineer.

MAJUG – Fall Meeting 2017

Thursday, September 28, 2017

8:30 am  |  Eastern Standard Time (New York, GMT-04:30)  |  4.0 hrs

[#2 · Accepted solution]

clear log() ;
/*****************************************************************
This JMP and R Integration script and example uses the

Alaska pipeline data consisting of in-field ultrasonic

measurements of the depths of defects in the Alaska pipeline.

The depth of the defects were then re-measured in the laboratory.
These measurements were performed in six different batches.
The data were analyzed to calibrate the bias of the field
measurements relative to the laboratory measurements.

In this analysis, the field measurement is the response variable
and the laboratory measurement is the predictor variable.

These data were provided by Harry Berger, who was at thetime a scientist for the Office of the Director of the Institute
of Materials Research (now the Materials Science and Engineering
Laboratory) of NIST.

These data were used for a study conducted for the Materials
Transportation Bureau of the U.S. Department of Transportation.

The source of the data came from the following link
http://www.itl.nist.gov/div898/handbook/pmd/section6/pmd621.htm

***************************************************************/

R Init() ;/* starts the connection to R */

R Submit( "\[
# load packages
library (lattice)

## Alaska pipeline ultrasonic calibration case study
## Create vector with dependent variable, field defect size
fdef = c(18,38,15,20,18,36,20,43,45,65,43,38,33,10,50,10,50,15,53,60,18, 38,15,
20,18,36,20,43,45,65,43,38,33,10,50,10,50,15,53,15,37,15, 18,11,35,20,40,50,36,
50,38,10,75,10,85,13,50,58,58,48,12,63,10, 63,13,28,35,63,13,45,9,20,18,35,20,
38,50,70,40,21,19,10,33,16,5, 32,23,30,45,33,25,12,53,36,5,63,43,25,73,45,52,
9,30,22,56,15,45)

## Create vector with independent variable, lab defect size
ldef = c(20.2,56.0,12.5,21.2,15.5,39.0,21.0,38.2,55.6,81.9,39.5,56.4,40.5,
14.3,81.5,13.7,81.5,20.5,56.0,80.7,20.0,56.5,12.1,19.6,15.5,38.8, 19.5,38.0,
55.0,80.0,38.5,55.8,38.8,12.5,80.4,12.7,80.9,20.5,55.0, 19.0,55.5,12.3,18.4,
11.5,38.0,18.5,38.0,55.3,38.7,54.5,38.0,12.0, 81.7,11.5,80.0,18.3,55.3,80.2,
80.7,55.8,15.0,81.0,12.0,81.4,12.5, 38.2,54.2,79.3,18.2,55.5,11.4,19.5,15.5,
37.5,19.5,37.5,55.5,80.0, 37.5,15.5,23.7,9.8,40.8,17.5,4.3,36.5,26.3,30.4,50.2,
30.1,25.5, 13.8,58.9,40.0,6.0,72.5,38.8,19.4,81.5,77.4,54.6,6.8,32.6,19.8,
58.8,12.9,49.0)

## Create vector with batch indicator
bat = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,
4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
5,6,6,6,6,6,6,6)

## Save data in a data frame and determine number of observations
Batch <- as.factor(bat)
df <- data.frame(fdef,ldef,Batch)
len <- length(Batch)
xax = "Lab Defect Size"
yax = "Field Defect Size"
title = "Alaska Pipeline Ultrasonic Calibration Data"
plot(ldef,fdef,xlab=xax,ylab=yax,main=title,col="blue")

## Generate conditional plot
trellis.device(new = TRUE, col = FALSE)
FIG = xyplot(fdef ~ ldef | Batch,
data=df,

main = title,
layout=c(3,2), col=4,
xlab=list(xax,cex=1.1), ylab=list(yax,cex=1.1),
strip=function(...) strip.default(...,strip.names=c(T,T)))
FIG

## Batch analysis
x = ldef
y = fdef
xydf = data.frame(x,y,Batch)
out = by(xydf,xydf$Batch,function(x) lm(y~x,data=x))
lapply(out,summary)

#> Residual standard error: 1.35 on 5 degrees of freedom
#> Multiple R-Squared: 0.9956, Adjusted R-squared: 0.9947
#> F-statistic: 1128 on 1 and 5 DF, p-value: 4.401e-07
## Save batch regression results
outs = sapply(out,summary)
outc = sapply(out,coef)
fitse = t(outs[6,])
fitse = c(fitse[[1]],fitse[[2]],fitse[[3]],fitse[[4]],fitse[[5]],fitse[[6]])
r2 = t(outs[8,])
r2 = c(r2[[1]],r2[[2]],r2[[3]],r2[[4]],r2[[5]],r2[[6]])
b0 = t(outc[1,])
b1 = t(outc[2,])
## Batch plots
par(mfrow=c(2,2))
id = c(1:length(b0))
xax2 = "Batch Number"
plot(id,r2,xlab=xax2,ylab="Correlation",ylim=c(.8,1), col="blue",pch=16,cex=1.25)
abline(h=mean(r2))

plot(id,b0[1,],xlab=xax2,ylab="Intercept",ylim=c(0,8), col="blue",pch=16,cex=1.25)
abline(h=mean(b0))

plot(id,b1[1,],xlab=xax2,ylab="Slope",ylim=c(.5,.9), col="blue",pch=16,cex=1.25)
abline(h=mean(b1))

plot(id,fitse,xlab=xax2,ylab="RESSD",ylim=c(0,7), col="blue",pch=16,cex=1.25)
abline(h=mean(fitse))
par(mfrow=c(1,1))
## Straight line regression analysis
out = lm(fdef~ldef)
summary(out)

## Residual 6-plot
par(mfrow=c(3,2))
plot(ldef,fdef,xlab="Lab Defect Size",
ylab="Field Defect Size",main="Field Defect Size vs Lab Defect Size")
abline(reg=out)
plot(ldef,out$residuals,ylab="Residuals",xlab="Lab Defect Size",
main="Residuals vs Lab Defect Size")
plot(out$fitted.values,out$residuals,ylab="Residuals",xlab="Predicted",
main="Residual vs Predicted")
plot(out$residuals[2:len],out$residuals[1:len-1],ylab="Residuals",
xlab="Lag 1 Residuals",main="Lag Plot")
hist(out$residuals,ylab="Frequency",xlab="Residuals",main="Histogram")
qqnorm(out$residuals,main="Normal Probability Plot")
## Generate plot of raw data with overlaid regression function
par(mfrow=c(1,1),cex=1.25)

plot(ldef,fdef,ylab="Field Defect Size",xlab="Lab Defect Size", col="blue")
abline(reg=out)
title("Alaska Pipeline Ultrasonic Calibration Data",line=2)
title("With Unweighted Line",line=1)

## Plot residuals versus lab defect size
par(mfrow=c(1,1),cex=1.25)
plot(ldef,out$residuals, xlab="Lab Defect Size", ylab="Residuals",
main="Alaska Pipeline Data Residuals - Unweighted Fit", cex=1.25, col="blue")
abline(h=0)

## Transformations of response variable
lnfdef = log(fdef)
sqrtfdef = sqrt(fdef)
invfdef = 1/fdef

## Plot transformed response variable
par(mfrow=c(2,2))
xax = "Lab Defect Size"
plot(ldef,fdef,xlab=xax,ylab="Field Defect Size",col="blue")
plot(ldef,sqrtfdef,xlab=xax,ylab="Sqrt(Field Defect Size)",col="blue")
plot(ldef,lnfdef,xlab=xax,ylab="ln(Field Defect Size)",col="blue")
plot(ldef,invfdef,xlab=xax,ylab="1/Field Defect Size",col="blue")
title(main="Transformations of Response Variable",outer=TRUE,line=-2)

## Transformations of predictor variable
lnldef = log(ldef)
sqrtldef = sqrt(ldef)
invldef = 1/ldef

## Plot transformed predictor variable
par(mfrow=c(2,2))
yax = "ln(Field Defect Size)"
plot(ldef,lnfdef,xlab="Lab Defect Size", ylab=yax, col="blue")
plot(sqrtldef,lnfdef,xlab="Sqrt(Lab Defect Size)", ylab=yax, col="blue")
plot(lnldef,lnfdef,xlab="ln(Lab Defect Size)", ylab=yax, col="blue")
plot(invldef,lnfdef,xlab="1/Lab Defect Size", ylab=yax, col="blue")
title(main="Transformations of Predictor Variable",outer=TRUE,line=-2)

]\" ) ;

/* Create JMP data table from the R data frames */
Rdf = R Get(df) ;
Rdf << NewDataView << Set Name("R Data frame Table") ;
/** Use Graph Builder on dtrdf to produce scatter plots **/
/*Set the Value ordering property for the Batch column so that it
matched displaying the Trellis Plot in R
*/
Rdf:Batch << Set Property( "Value Ordering", Eval ( {4, 5, 6, 1, 2, 3} ) ) <<
Value Labels(
{1 = "Batch: 1", 2 = "Batch: 2", 3 = "Batch: 3", 4 = "Batch: 4", 5 =
"Batch: 5", 6 = "Batch: 6"}
) <<
Use Value Labels( 1 )
;
/** use row states to change the plot markers from solid to blank to match R plot ***/
Rdf << Set Row States(
[128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128,
128, 128]
)
;

Rdf << Graph Builder(
Size( 457, 408 ),
Show Control Panel( 0 ),
Variables( X( :ldef ), Y( :fdef ), Wrap( :Batch, Levels per Row( 3 ) ) ),
Elements( Points( X, Y, Legend( 3 ) ) ),
SendToReport(
Dispatch(
{},
"graph title",
TextEditBox,
{Set Text( "Alsaka Pipeline Ultrasonic Calibration Data - JMP" )}
),
Dispatch( {}, "X title", TextEditBox, {Set Text( "Lab Defect Size " )} ),
Dispatch(
{},
"Y title",
TextEditBox,
{Set Text( "ln(Field Defect Size)" )}
),
Dispatch( {}, "Graph Builder", FrameBox( 1 ), {Marker Size( 1 )} ),
Dispatch( {}, "Graph Builder", FrameBox( 2 ), {Marker Size( 1 )} ),
Dispatch( {}, "Graph Builder", FrameBox( 3 ), {Marker Size( 1 )} ),
Dispatch( {}, "Graph Builder", FrameBox( 4 ), {Marker Size( 1 )} ),
Dispatch( {}, "Graph Builder", FrameBox( 5 ), {Marker Size( 1 )} ),
Dispatch( {}, "Graph Builder", FrameBox( 6 ), {Marker Size( 1 )} )
)
);

/* Produce plots and regression fits with JMP Analysis platforms */
/* Get the FIG graph from R and display in JMP
*/
R Submit( "\[
## Generate conditional plot
trellis.device(new = TRUE, col = FALSE)
FIG = xyplot(fdef ~ ldef | Batch,
data=df,

main = title,
layout=c(3,2), col=4,
xlab=list(xax,cex=1.1), ylab=list(yax,cex=1.1),
strip=function(...) strip.default(...,strip.names=c(T,T)))
FIG
]\" );

picture = R Get Graphics( "png" );
New Window("Alaska Pipeline Ultrasonic Calibration Data from R into JMP", picture) ;
Wait( 10 );

/* Box-Cox linearity plot */

R Submit( "\[
## Transformations of predictor variable
lnldef = log(ldef)
sqrtldef = sqrt(ldef)
invldef = 1/ldef

## Box-Cox linearity plot
for (i in (0:100)){
alpha = -2 + 4*i/100
if (alpha != 0){
tx = ((ldef**alpha) - 1)/alpha
temp = lm(lnfdef~tx)
temps = summary(temp)
if(i==0) {rsq = temps$r.squared
alp = alpha}
else {rsq = rbind(rsq,temps$r.squared)
alp = rbind(alp,alpha)}
}}
rcor = sqrt(rsq)
par(mfrow=c(1,1),cex=1.25)
plot(alp,rcor,type="l",xlab="Alpha",ylab="Correlation",
main="Box-Cox Linearity Plot ln(Field) Lab",
ylim=c(.6,1), col="blue")

#creaate Correlation and Alpha vectors into BCP data frame
Corr <- as.vector(rcor)
Alpha <- as.vector(alp)
BCP <- data.frame(Corr,Alpha)

]\" );
RBCP = R Get(BCP);
RBCP << NewDataView << Set Name("R Data Frame of Alpha-Correlation Table") ;

/* Graph Box-Cox Linearity Plot from Graph Builder */
RBCP << Graph Builder(
Size( 421, 362 ),
Show Control Panel( 0 ),
Variables( X( :Alpha ), Y( :Corr ) ),
Elements( Smoother( X, Y, Legend( 4 ) ) ),
SendToReport(
Dispatch(
{},
"graph title",
TextEditBox,
{Set Text( "Box-Cox Linearity Plot ln(Field) - JMP" )}
),
Dispatch( {}, "Y title", TextEditBox, {Set Text( "Correlation" )} )
)
);

BoxCox = R Get Graphics( "png" );
New Window("Box-Cox Linearity Plot ln(Field) Lab - From R", BoxCox) ;
Wait( 10 );

/*****

Rxydf <<
Graph Builder(
Size( 533, 454 ),
Show Control Panel( 0 ),
Variables( X( :x ), Y( :y ) ),
Elements( Points( X, Y, Legend( 3 ) ) ),
SendToReport(
Dispatch(
{},
"y",
ScaleBox,
{Label Row( Label Orientation( "Horizontal" ) )}
),
Dispatch(
{},
"graph title",
TextEditBox,
{Set Text( "Alaskan Pipeline Ultrasonic Calibration Data - JMP" )}
),
Dispatch( {}, "X title", TextEditBox, {Set Text( "Lab Defect Size" )} ),
Dispatch(
{},
"Y title",
TextEditBox,
{Rotate Text( "Left" ), Set Text( "Field Defect Size" )}
)
)
);
***/

/** Do a first initial regression straight-line model fit to the data

Rxydf << Bivariate(
Y( :y ),
X( :x ),
Fit Line( {Line Color( {57, 177, 67} )} ),
SendToReport(
Dispatch( {"Linear Fit"}, "Lack Of Fit", OutlineBox, {Close( 0 )} )
)
);

**/

R Term () ; /* terminates connection to R */

[#3]

Exploring JMP® Modeling Functionality Using Consumer Expenditure Data presentation and example data table with script. CE Public-Use Microdata accessible from: https://www.bls.gov/cex/pumd_data.htm  

Example JSL to process R BAS package: 

dt = Open(“C:\Desktop\Temp\JMP SESUG Materials\intrvw15\intrvw15\fmly2.jmp”);

R Init ( );

R Send (dt);

R Submit(“\[

library(BAS)

library(stats)

library(dplyr)

# Fit the BAS model

CUEXP_bma = bas.lm(log(TOTEXPCQ) ~ CUTENURE + VEHQ + log(HOUSCQ) + INCLASS + VEHQL + log(SHELTCQ) + TOTBATHQ + CHILDAGE + HIGH_EDU + CPIWX_FLAG_REF + CPIWX_FLAG_SPO + FAM_TYPE + FAM_SIZE + REF_RACE,

data = dt,

prior = "BIC",

modelprior = uniform()

)

CUEXP_bma

options(width = 100)

summary(CUEXP_bma)

　

# Root Mean Square Error

pred.HPM = predict(CUEXP_bma, newdata = dt, estimator = 'HPM', se.fit=TRUE)

# pred.HPM.rmse = sqrt(mean((exp(pred.HPM$fit) - dt$TOTEXPCQ)^2))

# cat('HPM_RMSE_TOTEXPCQ\n')

# pred.HPM.rmse

pred.HPM.rmse2 = sqrt(mean((pred.HPM$fit - log(dt$TOTEXPCQ))^2))

cat('HPM_RMSE_logTOTEXPCQ\n')

pred.HPM.rmse2

# Coverage Probability based on Confidence Interval

out = as.data.frame(cbind(exp(confint(pred.HPM)), TOTEXPCQ = dt$TOTEXPCQ ))

colnames(out)[1:2] = c('lwr', 'upr')

out %>% summarize(COVERAGE_CI = sum(TOTEXPCQ >= lwr & TOTEXPCQ <=upr) / n())

]\”);