Which test statistic should I use?

Question

I have collected data consisting of radiation dose from 2 different equipment. I have taken 3 readings (reps) on each exposure using both the equipment. The data looks like "at 5 mAs and 55 kVp" Equip1 readings in Gy (0.12, 0.11, 0.14) Equip2 readings (0.14, 0.15, 0.13). Like this I have range of mAs from 5 to 150 and also range of kVp from 55 to 120. I want to check if both the equipment are consistent in displaying the radiation dose. Any suggestions? Should I take mean of all 3 reps and then run a student's t test, or should I use split-split plot RCBD?

Answer 1

With respect to using a t-test, you would have to re-run it at each level of emitted radiation, and then find a good way to compare results across each level. Furthermore, with so few observations at each level, the results would likely be suspect on sample size issues, nevermind the general iffiness of trying to establish whether two groups are the same; you can certainly do it, as a non-statistically significant difference between them would indicate just that, but the traditional approach is to try to disprove the lack of a difference rather than prove its existence.

Running a 2-sided T test with the values you provided returns a p-value of 0.2, which would, in a very naive fashion, answer your question, of whether the two machines record the "different amount of radiation"; namely, assuming that we will only reject the null hypothesis at a significance level of 0.05 or less, then for this amount of emitted radiation, there are no differences in the recorded amounts. Again, massive, massive iffiness in the general approach here.

My advice? Take all your observations and put them through a linear model of recorded vs. actual, using a categorical variable to denote whether the reading was taken by machine 1 vs. machine 2. While you're at it, you might also want to consider interaction terms between radiation levels and machine identifier (e.g. machine 1 reading low-level radiation). Just don't get too carried away with it, lest the curse of dimensionality mess it all up for you.

So long as you perform your due diligence with respect to data prep and residual analysis doesn't throw up any red flags, no one should take issue with your results, whatever they might be.

See below for R commands and output of t-test on entire population:

Commands:

###Import csv datafile
data = data.table(read.csv('#:\\#####\\#####\\#####\\data.csv', header=T))
####re-arrange dataset to useful format for analysis
data[,9:10] = NULL
data11 = data[,c(1,2,3),with=F]
data12 = data[,c(1,2,4),with=F]
data13 = data[,c(1,2,5),with=F]
setnames(data11, 'R1', 'y')
setnames(data12, 'R2', 'y')
setnames(data13, 'R3', 'y')
data1 = rbind(data11, data12, data13)
data1$mch1 = 1

data21 = data[,c(1,2,6),with=F]
data22 = data[,c(1,2,7),with=F]
data23 = data[,c(1,2,8),with=F]
setnames(data21, 'R1.1', 'y')
setnames(data22, 'R2.1', 'y')
setnames(data23, 'R3.1', 'y')
data2 = rbind(data21, data22, data23)
data2$mch1 = 0

dataFinal = rbind(data1, data2)

###Population-level t-test
tTest = t.test(data1$y, data2$y)

Output:

Welch Two Sample t-test data: data1$y and data2$y t = -1.6514, df = 568.761, p-value = 0.09921 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.14518205 0.01255892 sample estimates: mean of x mean of y 0.4805459 0.5468575

Okay, so we can see that there might be a difference between the two machines, but we're not convinced just yet, that seemingly so-so p-value could be hiding a lot of stuff that got swept under the rug, so I ran a linear model to see what would shake out.

Long story short, there does indeed appear to be a difference between the two machines.

• The final model was built on a dataset that dropped the 15 (out of a total 588) observations with highest recorded readings as outliers (they had values higher than median+3*SD).

• The dependent variable was log transformed due to its exponential distribution (the transform wasn't brilliant, but I don't have the time to search for that one power transform that will do the trick), which required dropping an additional record with original reading of 0

• Machine 2 was used as the baseline

• Exponential terms of powers 2 and 3 were introduced for both continuous variables (mAs and kV)

• Interaction term between mAs and kV

Results, residual plots, and code:

Call:
lm(formula = log(y) ~ mAs + mAs2 + mAs3 + kV + kV2 + kV3 + mAs * 
kV + mAs3 * kV3 + mch1, data = dataFinal2[-1, ])
Residuals:     
Min       1Q   Median       3Q      Max 
-1.13542 -0.11617  0.01442  0.10731  2.43333 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.218e+01  9.480e-01 -12.852  < 2e-16***
mAs          1.105e-01  2.978e-03  37.120  < 2e-16***
mAs2        -1.136e-03  3.630e-05 -31.310  < 2e-16***
mAs3         3.891e-06  1.656e-07  23.500  < 2e-16***
kV           2.034e-01  3.425e-02   5.938 5.05e-09***
kV2         -1.523e-03  4.022e-04  -3.786  0.00017***
kV3          4.246e-06  1.534e-06   2.768  0.00582**
mch1        -5.763e-02  1.900e-02  -3.032  0.00254**
mAs:kV      -5.266e-05  2.329e-05  -2.261  0.02414*
mAs3:kV3    -1.157e-13  4.860e-14  -2.381  0.01761*
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2272 on 562 degrees of freedom
Multiple R-squared:  0.9657,    Adjusted R-squared:  0.9651 
F-statistic:  1757 on 9 and 562 DF,  p-value: < 2.2e-16

###Remove outliers and take log transform
dataFinal2 = dataFinal[(dataFinal$y<median(dataFinal$y)+3*sd(dataFinal$y))
dataFinal2$mAs2 = dataFinal2$mAs**2
dataFinal2$mAs3 = dataFinal2$mAs**3
dataFinal2$kV2 = dataFinal2$kV**2
dataFinal2$kV3 = dataFinal2$kV**3
###Build model and view results
lm2 = lm(log(y) ~ mAs+mAs2+mAs3+kV+kV2+kV3+mAs*kV+mAs3*kV3+mch1, dataFinal2[-1,])
summary(lm2)
plot(lm2)

That's where I'm leaving it. The model and its conclusions can still be improved, such as through the use of a better transformation for the dependent variable, the iterative removal of influential and high-leverage observations, to say nothing of testing its reliability thru the use of a test set (or even cross-validation, if you're really gung-ho), but you'll have to stop at some point and this is mine. Hope this helps.