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 */
... View more