Electronic Supplementary Material
Accreditation and Quality Assurance, 2012, Springer, Heidelberg
Experimental design and data evaluation considerations
for comparisons of reference materials
David L. Duewer, Hugo Gasca Aragon, Katrice A. Lippa, Blaza Toman
National Institute of Standards and Technology (NIST)
Gaithersburg, MD 20899 USA
Disclaimer
Certain commercial software is identified in this report to specify the experimental procedure as completely as possible. In no case does such identification imply a recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the software is necessarily the best available for the purpose.
Contents
ANOVA-based tests for determining degrees of freedom...... S-
Software for estimation of expanded uncertainties...... S-
Excel “ANOVA: Single Factor”...... S-
SAS MIXED procedure...... S-
One-factor with replication...... S-
Two-factor nested...... S-
R lme4...... S-
One-factor with replication...... S-
Two-factor nested...... S-
OpenBUGS...... S-
One-Factor with replication...... S-
Two-factor nested...... S-
Exemplar GDR analysis...... S-
FREML...... S-
Data and commands...... S-
Tabular results...... S-
Graphical results...... S-
XLGENLINE...... S-
Data, commands, and tabular results...... S-
Graphical results...... S-
RegViz...... S-
Data and commands...... S-
Tabular results...... S-
Graphical results...... S-
References...... S-
RM Comparison Studies (Electronic Supplementary Material)
ANOVA-based tests for determining degrees of freedom
For a one-factor design with replication, the estimated number of degrees of freedom are
and for a two-factor nested design they are
For these designs, assuming random effects, there are two frequentist tests for the presence of a variance component based on the mean square errors: the classical significance test [[1]] and the less common test for a variance component estimate to be positive [[2]]. The 95 % level of confidence test for a variance component classical estimate uses the criteria:
where Mb is the between-group mean square, Mw is the within-group mean square, and Mr is the residual mean square. This test requires the variance component estimate being large enough for statistical significance. It will thus confirm the obvious outcome for small samples with large random effects but is insensitive to relatively small variance components.
The test for a variance component estimate to be positive uses the criteria:
.
This addresses whether the random effects are likely to be present regardless of the magnitude of the variance component estimate. This appears to be the test used in the SAS MIXED procedure. See below for an implementation of this test in R.
Software for estimation of expanded uncertainties
The following data and analyses procedures are intended as examples to guide use of the software systems and/or to benchmark analyses using other software systems. While many values in the different sets of results are useful for estimating variance components, the critical values for estimating expanded uncertainties on a mean result are in red font. The information provided by the Excel, SAS, and R implementations of the frequentist estimates discussed below all result in the same estimated values, although sometimes requiring additional calculations.
The same data are used in all examples. Data for two materials are provided for the one-factor with replication design, with and without appreciable between-campaign variance. Data for four exemplar materials are provided for the two-factor nested design, with and without appreciable between-aliquot and/or between-campaign variance.
Excel “ANOVA: Single Factor”
The Microsoft Office Excel 2003 “ANOVA: Single Factor” is an extremely easy-to-use tool that is adequate for evaluating one-factor designs with replication, but does require additional effort to interpret. The tool is provided in the “Analysis ToolPak” Excel add-in [[3]], an optional suite of data analysis tools freely available for most if not all versions of Excel. Most basic data analysis systems support one-factor ANOVA and it’s readily doable “by hand” [[4]].
Input data and exemplar command
Material A: No appreciable between-campaign variance
ANOVA: Single FactorSUMMARY
Groups / Count / Sum / Average / Variance
Cmp1 / 3 / 152584 / 50861.3 / 36601.3
Cmp2 / 3 / 152842 / 50947.3 / 1192.3
ANOVA
Source of Variation / SS / df / MS / F / P-value / F crit
Between Groups / 11094.0 / 1 / 11094.0 / 0.587 / 0.4863 / 7.709
Within Groups / 75587.3 / 4 / 18896.8
Total / 86681.3 / 5
The “P-value” of 0.49 indicates that the “Between-Groups” (here, between-campaign) contribution to the variance is not significant. The standard uncertainty of the mean is therefore best calculated as the standard deviation/square root of the number of measurements. This tool provides neither the mean nor the SD, but they are easily calculated using Excel’s “AVERAGE(.)” and “STDEV(.)” functions. The mean is 50904.3, the standard deviation is 131.7, the standard uncertainty of the mean is , and the ±U95 is therefore t0.05,5 × 54 = 2.6 × 54 = 138.
Results for Material B: appreciable between-campaign variance
ANOVA: Single FactorSUMMARY
Groups / Count / Sum / Average / Variance
Cmp1 / 3 / 70763 / 23587.7 / 158.3
Cmp2 / 3 / 70242 / 23414.0 / 3079.0
Cmp3 / 3 / 70942 / 23647.3 / 297.3
ANOVA
Source of Variation / SS / df / MS / F / P-value / F crit
Between Groups / 88164.7 / 2 / 44082.3 / 37.414 / 0.0004 / 5.143
Within Groups / 7069.3 / 6 / 1178.2
Total / 95234.0 / 8
The “P-value” of 0.0004 indicates that there is a significant between-campaign contribution to the variance. The standard uncertainty on the mean is and v = 2. The ±U95 is thus estimated to be t0.05,2 × 70 = 4.3 × 70 = 301.
SAS MIXED Procedure
SAS [[5]] is an extensively documented and validated commercial data analysis system. The MIXED procedure [[6]] provides analysis of linear mixed-effects models.
It requires a good deal of experience and expertise to efficiently set up a SAS analysis and to select among the many options that are available for most of its component procedures, MIXED included. SAS is very definitely not freeware, but many institutions and universities do have licensed access.
One-factor with replication
Using SAS to analyze a one-factor completely balanced design is rather like using a sledgehammer to kill ants, but it does provide a “definitive” benchmark for less comprehensive (and less expensive) analysis systems.
Data and commands for the MIXED procedure
data CCQM_1F_Exemplar;
input material $char14. cmp rep y;
datalines;
MATERIAL A 1 1 50742
MATERIAL A 1 2 50760
MATERIAL A 1 3 51082
MATERIAL A 2 1 50912
MATERIAL A 2 2 50949
MATERIAL A 2 3 50981
MATERIAL B 1 1 23602
MATERIAL B 1 2 23576
MATERIAL B 1 3 23586
MATERIAL B 2 1 23431
MATERIAL B 2 2 23353
MATERIAL B 2 3 23459
MATERIAL B 3 1 23667
MATERIAL B 3 2 23645
MATERIAL B 3 3 23632
;
proc sort data= CCQM_1F_Exemplar;
by material;
run;
proc mixed data=CCQM_1F_Exemplar asycov noitprint;
class cmp rep;
model y = / solution ddfm=sat;
random cmp;
by material;
run;
Material A, no appreciable between-campaign variance
------material=MATERIAL A ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_1F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 2 1 2
rep 3 1 2 3
Dimensions
Covariance Parameters 2
Columns in X 1
Columns in Z 2
Subjects 1
Max Obs Per Subject 6
Number of Observations
Number of Observations Read 6
Number of Observations Used 6
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 0
Residual 17336
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2
1 cmp
2 Residual 1.2022E8
Fit Statistics
-2 Res Log Likelihood 64.8
AIC (smaller is better) 66.8
AICC (smaller is better) 68.1
BIC (smaller is better) 65.5
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 50904 53.7529 5 947.01 <.0001
The “Intercept” (mean) is 50904 with a “Standard Error” (standard uncertainty) of 54 and v = 5. The ±U95 estimate for the mean is thus t0.05,5 × 54 = 2.6 × 54 = 138.
Material B, appreciable between-campaign variance
------material=MATERIAL B ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_1F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 3 1 2 3
rep 3 1 2 3
Dimensions
Covariance Parameters 2
Columns in X 1
Columns in Z 3
Subjects 1
Max Obs Per Subject 9
Number of Observations
Number of Observations Read 9
Number of Observations Used 9
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 14338
Residual 1167.44
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2
1 cmp 2.1693E8 -151436
2 Residual -151436 454309
Fit Statistics
-2 Res Log Likelihood 88.7
AIC (smaller is better) 92.7
AICC (smaller is better) 95.1
BIC (smaller is better) 90.9
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 23550 70.0636 2 336.12 <.0001
The mean is 23550 with a standard uncertainty of 70 and v = 2. The ±U95 estimate for the mean is thus t0.05,2 × 70 = 4.3 × 70 = 301.
Two-factor nested
Data and commands for the MIXED procedure
data CCQM_2F_Exemplar;
input material $char14. campaign aliquot rep y;
datalines;
MATERIAL C 1 1 1 37.451
MATERIAL C 1 1 2 38.416
MATERIAL C 1 2 1 37.720
MATERIAL C 1 2 2 38.268
MATERIAL C 2 1 1 38.054
MATERIAL C 2 1 2 37.386
MATERIAL C 2 2 1 37.934
MATERIAL C 2 2 2 38.484
MATERIAL D 1 1 1 34.045
MATERIAL D 1 1 2 34.562
MATERIAL D 1 2 1 34.031
MATERIAL D 1 2 2 33.249
MATERIAL D 2 1 1 34.421
MATERIAL D 2 1 2 34.320
MATERIAL D 2 2 1 33.852
MATERIAL D 2 2 2 34.512
MATERIAL E 1 1 1 27.370
MATERIAL E 1 1 2 27.436
MATERIAL E 1 2 1 27.170
MATERIAL E 1 2 2 27.467
MATERIAL E 2 1 1 27.508
MATERIAL E 2 1 2 27.429
MATERIAL E 2 2 1 27.431
MATERIAL E 2 2 2 27.651
MATERIAL F 1 1 1 7.105
MATERIAL F 1 1 2 6.983
MATERIAL F 1 2 1 7.067
MATERIAL F 1 2 2 7.042
MATERIAL F 2 1 1 7.132
MATERIAL F 2 1 2 7.096
MATERIAL F 2 2 1 7.172
MATERIAL F 2 2 2 7.422
;
proc sort data=CCQM_2F_Exemplar;
by material;
run;
proc mixed data=CCQM_2F_Exemplar asycov noitprint;
class campaign aliquot rep;
model y = / solution ddfm=sat;
random campaign aliquot(campaign);
by material;
run;
Material C, no appreciable between-aliquot or -campaign variance
------material=MATERIAL C ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_2F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 2 1 2
alq 2 1 2
rep 2 1 2
Dimensions
Covariance Parameters 3
Columns in X 1
Columns in Z 6
Subjects 1
Max Obs Per Subject 8
Number of Observations
Number of Observations Read 8
Number of Observations Used 8
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 0
alq(cmp) 0
Residual 0.1761
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2 CovP3
1 cmp
2 alq(cmp)
3 Residual -0.00274 0.008863
Fit Statistics
-2 Res Log Likelihood 9.8
AIC (smaller is better) 11.8
AICC (smaller is better) 12.6
BIC (smaller is better) 10.5
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 37.9641 0.1484 7 255.86 <.0001
The mean is 37.96 with a standard uncertainty of 0.15 and v = 7. The ±U95 is thus estimated to be t0.05,7 × 0.15 = 2.4 × 0.15 = 0.36. Note that the between-campaign and between-aliquot covariance estimates are both zero.
Material D, appreciable between-aliquot variance
------material=MATERIAL D ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_2F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 2 1 2
alq 2 1 2
rep 2 1 2
Dimensions
Covariance Parameters 3
Columns in X 1
Columns in Z 6
Subjects 1
Max Obs Per Subject 8
Number of Observations
Number of Observations Read 8
Number of Observations Used 8
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 0
alq(cmp) 0.02741
Residual 0.1656
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2 CovP3
1 cmp
2 alq(cmp) 0.01152 -0.00685
3 Residual -0.00685 0.01371
Fit Statistics
-2 Res Log Likelihood 10.2
AIC (smaller is better) 14.2
AICC (smaller is better) 17.2
BIC (smaller is better) 11.6
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 34.1240 0.1660 3 205.59 <.0001
The mean is 34.12 with a standard uncertainty of 0.17 and v = 3. The ±U95 is thus estimated to be t0.05,3 × 0.17 = 3.2 × 0.17 = 0.54. Note that the between-campaign covariance estimate is zero but that the between-aliquot estimate is greater than zero.
Material E, appreciable between-campaign variance
------material=MATERIAL E ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_2F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 2 1 2
alq 2 1 2
rep 2 1 2
Dimensions
Covariance Parameters 3
Columns in X 1
Columns in Z 6
Subjects 1
Max Obs Per Subject 8
Number of Observations
Number of Observations Read 8
Number of Observations Used 8
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 0.006784
alq(cmp) 0
Residual 0.01433
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2 CovP3
1 cmp 0.000219 -0.00002
2 alq(cmp)
3 Residual -0.00002 0.000068
Fit Statistics
-2 Res Log Likelihood -6.7
AIC (smaller is better) -2.7
AICC (smaller is better) 0.3
BIC (smaller is better) -5.3
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 27.4328 0.07200 1 381.03 0.0017
The mean is 27.433 with a standard uncertainty of 0.072 and v = 1. The ±U95 is thus estimated to be t0.05,1 × 0.072 = 12.7 × 0.072 = 0.91. Note that the between-aliquot covariance estimate is 0 but the between-campaign estimate is greater than 0.
Material F, appreciable between-campaign and -aliquot variance
------material=MATERIAL F ------
The Mixed Procedure
Model Information
Data Set WORK.CCQM_2F_EXEMPLAR
Dependent Variable y
Covariance Structure Variance Components
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
cmp 2 1 2
alq 2 1 2
rep 2 1 2
Dimensions
Covariance Parameters 3
Columns in X 1
Columns in Z 6
Subjects 1
Max Obs Per Subject 8
Number of Observations
Number of Observations Read 8
Number of Observations Used 8
Number of Observations Not Used 0
Covariance Parameter
Estimates
Cov Parm Estimate
cmp 0.008007
alq(cmp) 0.003443
Residual 0.009913
Asymptotic Covariance Matrix of Estimates
Row Cov Parm CovP1 CovP2 CovP3
1 cmp 0.000316 -0.00004
2 alq(cmp) -0.00004 0.000083 -0.00002
3 Residual -0.00002 0.000049
Fit Statistics
-2 Res Log Likelihood -7.7
AIC (smaller is better) -1.7
AICC (smaller is better) 6.3
BIC (smaller is better) -5.6
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 7.1274 0.07812 1 91.23 0.0070
The mean is 7.127 with a standard uncertainty of 0.078 and v = 1. The ±U95 is thus estimated to be t0.05,1 × 0.078 = 12.7 × 0.078 = 0.99. Note that both the between-aliquot and between-campaign covariance estimates are greater than zero.
R lme4
R is an extremely powerful open-source environment for statistical computing and graphics [[7]]. The lme4 package [[8]] provides analysis of linear mixed-effects models.
Successfully getting R to run on your computer requires some “computer smarts” and successfully using R requires some programming and a bit of patience. However, R has become the tool of choice for many number-crunching tasks and learning to use it effectively is well worth the effort for anyone seriously interested in the analysis of data.
One-factor with replication
Input data and exemplar command
#======
# One Factor With Replication ANOVA for CCQM_1F_Exemplar Data #
#======
library(lme4)
# CCQM_1F_Exemplar Data
n<- 2
n.cmp<- 3
n.rep<- 3
camp<- as.factor(c(1,1,1,2,2,2,3,3,3))
d <- matrix(c(
50742, 50760, 51082, 50912, 50949, 50981, NA, NA, NA,
23602, 23576, 23586, 23431, 23353, 23459, 23667, 23645, 23632), n, n.cmp*n.rep, byrow=T)
material<- c("MATERIAL A","MATERIAL B")
value<- rep(0, n)
u.value<- rep(0, n)
u.camp<- rep(0, n)
u.replicate<- rep(0, n)
dof<- rep(0,n)
for (j in 1:n) {
y_ij <- d[j, c(1:(n.cmp*n.rep))]
value[j]<- mean(y_ij, na.rm=T)
dat <- data.frame(x=y_ij, camp=camp)
lmer.res<- lmer(x~1+(1|camp), data=dat, na.action=na.omit)
sigmas<- as.double(summary(lmer.res)@REmat[,4])
u.camp[j]<- sigmas[1]
u.replicate[j]<- sigmas[2]
u.value[j]<- sqrt(vcov(lmer.res)[1,1])
if (u.camp[j]==0) {
dof[j]<- length(y_ij[!is.na(y_ij)])-1
} else {
dof[j]<- length(levels(as.factor(camp)))-1
}
print( sprintf( "%14s %7.4f %7.4f %7.4f %7.4f %7.4f %2.0f",
material[j], value[j], u.value[j], qt(0.975,dof[j])*u.value[j], u.camp[j], u.replicate[j], dof[j] ) )
}
Results
[1] " MATERIAL A 50904.3333 53.7541 138.1792 0.0000 131.6700 5"
[1] " MATERIAL B 23550.1111 70.0638 301.4603 119.7400 34.1680 2"
The ±U95 estimates are 138 and 301.
Two-factor nested
Input data and exemplar command
#======
# Two-Factor Nested ANOVA for CCQM_2F_Exemplar Data #
======
library(lme4)
# CCQM_2F_Exemplar Data
n<- 4
n.cmp<- 2
n.ali<- 2
n.rep<- 2
camp<- as.factor(c(1,1,1,1,2,2,2,2))
aliquot<- as.factor(c(1,1,2,2,3,3,4,4))
d <- matrix(c(
37.451, 38.416, 37.720, 38.268, 38.054, 37.386, 37.934, 38.484,
34.045, 34.562, 34.031, 33.249, 34.421, 34.320, 33.852, 34.512,
27.370, 27.436, 27.170, 27.467, 27.508, 27.429, 27.431, 27.651,
7.105, 6.983, 7.067, 7.042, 7.132, 7.096, 7.172, 7.422
), n, n.cmp*n.ali*n.rep, byrow=T)
material<- c("MATERIAL C","MATERIAL D","MATERIAL E","MATERIAL F")
value<- rep(0, n)
u.value<- rep(0, n)
u.camp<- rep(0, n)
u.aliquot<- rep(0, n)
u.replicate<- rep(0, n)
dof<- rep(0,n)
for (i in 1:n) {
y_ijk <- d[i, 1:(n.cmp*n.ali*n.rep)]
value[i]<- mean(y_ijk)
dat <- data.frame(x=y_ijk, camp=camp, aliquot=aliquot)
# get the MS terms to test for the presence of variance components
tmp.aov<- summary(aov(x~camp*aliquot, data=dat))
MS<- tmp.aov[[1]][,3]
test.negative.var.camp<- pf(MS[2]/MS[1], n.cmp-1, n.cmp*(n.ali-1))>0.5
test.negative.var.aliquot<- pf(MS[3]/MS[2], n.cmp*(n.ali-1), n.cmp*n.ali*(n.rep-1))>0.5
lmer.res<- lmer(x~1+(1|camp)+(1|aliquot), data=dat)
sigmas<- as.double(summary(lmer.res)@REmat[,4])
u.camp[i]<- sigmas[2]
u.aliquot[i]<- sigmas[1]
u.replicate[i]<- sigmas[3]
u.value[i]<- sqrt(vcov(lmer.res)[1,1])
if (test.negative.var.camp) {
if (test.negative.var.aliquot) {
dof[i]<- n.cmp*n.ali*n.rep-1
} else {
dof[i]<- n.cmp*n.ali-1
}
} else {
dof[i]<- n.cmp-1
}
print( sprintf( "%14s %7.4f %7.4f %7.4f %2.0f %7.4f %7.4f %7.4f",
material[i], value[i], u.value[i], qt(0.975,dof[i])*u.value[i], dof[i], u.camp[i], u.aliquot[i],
u.replicate[i] ) )
}
Results
[1] " MATERIAL C 37.9641 0.1484 0.3509 7 0.0000 0.0000 0.4197"
[1] " MATERIAL D 34.1240 0.1660 0.5283 3 0.0000 0.1656 0.4069"
[1] " MATERIAL E 27.4327 0.0720 0.9148 1 0.0824 0.0000 0.1197"
[1] " MATERIAL F 7.1274 0.0781 0.9927 1 0.0895 0.0587 0.0996"
The ±U95 estimates are 0.35, 0.53, 0.91, and 0.99.
OpenBUGS
OpenBUGS is an empirical Bayesian analysis freeware system that runs in several environments including R.. OpenBUGS is an actively developed and maintained lineal descendant of the pioneering (but now static) WinBUGS system [[9]]. While we have occasionally encountered problems running this code in WinBUGS, we have not observed problems with OpenBUGS.
One-Factor with replication
Model
### 2-level Bounded Nested ANOVA for CCQM_1F_Exemplar###
# ******** Specified parameters *******
# n0 number of materials
# n1 number campaigns per material
# n2 number repeat measurements per campaign
# RCRM[n0] vector of u(CRM)/crm
# Ufact....bounding factor for RCRMs
# y[n0,n1,n2] matrix of measurements
#
# ******** Interesting Estimated parameters *******
# m0[n0] vector of material-response means
# r1 scalar instrumental imprecision relative SD
# s0[n0] vector of between-campaign SDs
#
# ******** Intermediate estimated parameters *******
# m1[n0,n1] vector of campaign-response means
# sbnd[n0] vector of prior distribution on s0: RCRM*m0
# t0[n0] vector of inverse between-variances: 1/s0^2
# t1[n0] vector of inverse within-variances: 1/(r1*m0)^2
#
# ******** Priors *******
# m0 Normal, very broad centered at 0
# m1 Normal, data defined
# r1 Uniform, between 0 and 10% relative
# s0 Uniform, between 0 and RM’s assigned uncertainty
# y Normal, data defined
#
# The only tricky bits here are:
# 1) “sbnd” is used to limit the variability of the response to be
# no greater than that of the assigned uncertainty. They’re
# computed from the relative assigned uncertainties – see RCRM.
# 2) The relative instrumental imprecision is assumed to be the same
# for all materials
#
Model{
r1 ~ dunif(0.0,0.1)
for( i0 in 1 : n0) {
m0[i0] ~ dnorm(0,0.0000001)
sbnd[i0] <- Ufact*RCRM[i0]*abs(m0[i0])
s0[i0] ~ dunif(0.0,sbnd[i0])
t1[i0] <- 1/pow(r1*m0[i0],2)
t0[i0] <- 1/pow(s0[i0],2)
for( i1 in 1 : n1 ) {m1[i0,i1] ~ dnorm(m0[i0],t0[i0])
for( i2 in 1 : n2 ) {y[i0,i1,i2] ~ dnorm(m1[i0,i1],t1[i0])}}}}
Data and initialization parameters
# CCQM_1F_Exemplar data
list(n0=2,n1=3,n2=3,Ufact=1,
RCRM=c(0.00048,0.0058),
y=structure(.Data=c(50742,50760,51082,50912,50949,50981,NA,NA,NA,
23601,23576,23586,23431,23352,23459,23666,23644,23632),.Dim = c(2,3,3)))
# CCQM_1F_Exemplar Initialization values (not required but speeds things up)
list(m0=c(5.1E+04,2.4E+04))
Results
Mean sd MC_errorval2.5pcmedian val97.5pcstartsample
m0[1]50889.052.6291.160950780.050890.050991.010001100000
m0[2]23541.061.7550.2105223415.023542.023662.010001100000
The ±U95 for material “A” (here labeled m0[1]) is estimated to be (50991-50780)/2 = 106 and for material “B” (m0[2]) is (23662-23415)/2 = 124, both shorter than the 139 and 301 of the frequentist evaluations.
Two-factor nested
Model
### 2-Factor Nested ANOVA for CRM KCs###
# ******** Specified parameters *******
# n0 number of materials
# n1 number campaigns per material
# n2 number aliquots per campaign
# n3 number replicates per aliquot
# RCRM[n0] vector of u(CRM)/crm
# Ufact bounding factor for RCRMs
# s1hig highest allowable value for relative between-aliquot SD