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 */