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?


1 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:


###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)


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:

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

          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

enter image description here enter image description here enter image description here enter image description here

###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,])

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.

  • $\begingroup$ Yes Welch Two Sample t-test worked for me $\endgroup$
    – Haider Ali
    Apr 16, 2015 at 15:45
  • $\begingroup$ @Haider Ali - Glad to hear it. I would still recommend using a linear model for your analysis though. Less grunt work required and a cleaner conclusion at the end. $\endgroup$
    – habu
    Apr 16, 2015 at 15:49
  • $\begingroup$ I am not getting the desired results using linear model. Can you help me out? Should I send you the data in excel file? $\endgroup$
    – Haider Ali
    Apr 16, 2015 at 16:23
  • $\begingroup$ Why not? Send it over and i'll have a look. $\endgroup$
    – habu
    Apr 16, 2015 at 16:56
  • $\begingroup$ how do I send it to you? Email? I don't see any attachment option $\endgroup$
    – Haider Ali
    Apr 16, 2015 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.