Question

1. Researchers at National Institutes of Standards and Technology (NIST) collected pipeline data on ultrasonic measurements of the depth of defects in the Alaska pipeline in the field. The depth of the defects were then remeasured in the lab- oratory. These measurements were performed in six different batches. It turns out that this batch effect is not significant and so can be ignored in the analysis that follows. The laboratory measurements are more accurate than the in-field measurements, but more time consuming and expensive. We want to develop a regression equation for correcting the in-field measurements. (a) Fit a regression model Lab Field. Check for non-constant variance. (b) We wish to use weights to account for the non-constant variance. Here we split the range of Field into 12 groups of size nine (except for the last group which has only eight values). Within each group, we compute the variance of Lab as varlab and the mean of Field as meanfield. Supposing pipeline is the name of your data frame, the following R code will make the needed computations: > i <- order (pipeline $Field) > npipe <- pipeline [i,] > ff <- gl (12,9) (-108] > meanfield <- unlist (lapply (split (npipe$Field, ff), mean)) > varlab unlist (lapply (split (npipe$Lab, ff),var)) <-

Answer 1

(a) Fit a regression model Lab ~ Field. Check for non-constant variance.

> rm(list = ls())
> require(nlme)
> data(pipeline, package="faraway")
> lm.fit <- lm(Lab ~Field, data = pipeline)
> summary(lm.fit)

> plot(pipeline$Field,pipeline$Lab, col=as.factor(pipeline$Batch),main = TeX("$Lab \\sim Field$"),pch='*', cex=1.5)
> legend("topright",title="Batch",legend= unique(pipeline$Batch), fill=1:length(pipeline$Batch) )
> abline(coef(lm.fit),lty=5)
> plot(residuals(lm.fit) ~ Field, pipeline, main ="Residuals versus log(time) for simple linear model")

(b) We wish to use weights to account for the non-constant variance.Here we split the range of Field into 12 groups of size nine (except for the last group which has only eight values). Within each group, we compute the variance of Lab as varlab and the mean of Field as meanfield. Supposing pipeline is the name of your data frame, the following R code will make the needed computations:

> i <- order(pipeline$Field)
> npipe <- pipeline[i,]
> ff <- gl(12,9)[-108]
> meanfield <- unlist(lapply(split(npipe$Field,ff),mean))
> varlab <- unlist(lapply(split(npipe$Lab,ff),var))
Suppose we guess that the the variance in the response is linked to the predictor in the following way:

var(Lab) = a0 Fielda1

Regress log(varlab) on log(meanfield) to estimate a0 and a1. (You might choose to remove the last point.) Use this to determine appropriate weights in a WLS fit of Lab on Field. Show the regression summary.

> plot(log10(varlab) ,log10(meanfield))
> df<-data.frame(cbind(as.numeric(meanfield),as.numeric(varlab)))
> df <- df[-c(12),]
> lm.var.model <- lm(log(varlab) ~ log(meanfield) , data= df)
> summary(lm.var.model)

As $var(Lab) = a0 \;Field^{a1}$ we have that $log(var(Lab)) = log(a0) + a1 \; log(Field)$

calculating weight vector with variances obtained from our model.

> pipeline<- pipeline[with(pipeline, order(Field)), ]
> a0 <-10^summary(lm.var.model)$coefficients[1]
> a1 <-summary(lm.var.model)$coefficients[2]
> var.lab <- a0 * pipeline$Field^a1
> se.lab <- sqrt(var.lab)

Ploting data with error bars from varience

> plot(pipeline$Field,pipeline$Lab, pch='*', ylim = c(0,100), main = "Field versus Lab with error bars")
> arrows(pipeline$Field, pipeline$Lab-se.lab,pipeline$Field, pipeline$Lab+se.lab, length=0.05, angle=90, code=3)

> pipeline$lab.var <- var.lab
> wlm.fit <- gls( Lab ~ Field, data=pipeline, weights = ~ var.lab)
> summary(wlm.fit)
> plot(pipeline$Field,pipeline$Lab, main = TeX("$Lab \\sim Field$ weighted regression with var model as weight"))
> abline(coef(wlm.fit),lty=5)
> plot(residuals(wlm.fit) ~ Field,pipeline, main ="Residuals versus Field for weighted regression ")

(c) An alternative to weighting is transformation. Find transformations on Lab and/or Field so that in the transformed scale the relationship is approximately linear with constant variance. You may restrict your choice of transformation to square root, log and inverse.

> lm.fit <- lm((Lab)^0.5 ~log(Field), data = pipeline)
> summary(lm.fit)

> plot(log(pipeline$Field),pipeline$Lab^0.5, col=as.factor(pipeline$Batch),main = TeX("$Lab \\sim Field$"),pch='*', > > cex=1.5)
> legend("topright",title="Batch",legend= unique(pipeline$Batch), fill=1:length(pipeline$Batch) )
> abline(coef(lm.fit),lty=5)
> plot(residuals(lm.fit) ~ Field, pipeline, main ="Residuals versus log(time) for simple linear model")