Developing and Communicating Effective Process Control Algorithms for the Clinical Laboratory

George V. Woodrow III

Background

The modern clinical lab is a highly automated enterprise. Typical instruments can process from 200 to 1000 samples per hour for 1 to 50 assays per sample.

QC has hardly changed from the days when techs could do 20 assays an hour manually.

The measurement is the ‘product’.

Quality control is done using Controls

Controls are surrogates for samples, mimicing their performance, but with known analyte concentration.

Controls are typically run at a rate of 2% to 4%, but can be run just once a day.

Assays have Total Allowable Error Requirements ((Subscript[TE, a])), set by the government or professional licensing agencies.

The Model

Two properties affect the accuracy of the measurement: bias and imprecision. Both typically vary with time.

Bias is the difference between the mean of a large number of controls and target value.

bias = μ - target

Total Error (TE) is modeled by the formula:

TE = |bias| + (Subscript[σ, m])×SD

where (Subscript[σ, m])is the sigma multiplier, and determines the maximum allowable defect rate.

defect Rate = 1 - (∫_(-∞)^(Subscript[σ, m]))(1/(Sqrt[2π]))(e^((-(x^2))/2))dx

This is just the area under the tail of the Gaussian distribution from (Subscript[σ, m])to ∞. Here is the sigma multiplier for various defect rates.

Grid[{{"Defect Rate %","(Subscript[σ, m])"},{5.0,1.65},{2.0,2.06},{1.0,2.33},{0.5,2.58},{0.1,3.1}}, Spacings -> {3,1}, Alignment -> {{Left, Center}, Automatic}, Frame -> All, Background -> {None, {{{LightGray, Gray}}, {1 -> LightYellow}}}, BaseStyle -> "Text"]

If TE ≤ (Subscript[TE, a]), the QC is acceptable. Otherwise, there is a QC failure and the samples must be reprocessed.

Pass/Fail Quality Control

Traditionally, there are not enough controls to compute the mean and SD of the controls accurately. Instead control limits are established and rules set that determine if the QC is acceptable. The most common rule is that with 3 controls, the process fails if one control is beyond the limits. Sometimes, rules are set so that a failue occurs if 2 of 3 controls are beyond limits.

Limits are set empirically. Sometimes, the limits are a simple multiple of the monthly control SD, and the limits are expressed by rules: 1-2s, or 1-3s, meaning that the process will fail if one control is beyond a limit computed by multiplying the monthly SD by 2, or by 3 in the second case. More typically, limits are modified empirically, making them tighter in response to client complaints or proficiency testing failures, and making them wider to lower the rate of reprocessing.

The process SD is assumed to be constant for all these rules.

Analytic Determination of Optimal QC Acceptance Limits

The ideal QC function has a probability of error detection of 0 when TE ≤ (Subscript[TE, a]), and 1 then TE > (Subscript[TE, a]). If we assume that the process SD is either constant, or changes very slowly, we can determine the critical bias for the assay and control. The critical bias is the bias when the Total Error equals the Allowable Total Error:

(Subscript[TE, a]) = TE = |CriticalBias| + (Subscript[σ, m])×SD
CriticalBias = (Subscript[TE, a]) - (Subscript[σ, m])×SD

Note that the larger the SD, the smaller the Critical Bias. The ideal error detection curve can be goven in terms of the critical bias:

criticalBias[tea_, sm_, procSD_] := tea - sm *procSD

pEDIdeal[x_, critBias_] := (Piecewise[{{ 0 , Abs[x] <= critBias}, {1, Abs[x] > critBias}}])

Even though the bias and SD of the process is unknown, the process does have these parameters. We can determine the probability that a single control is with a set of limits for a given bias and SD. It is just the area under the Gaussian distribution with parameters bias = μ - target and σ between ±limit:

(Subscript[P, s])=(1/(σ (Sqrt[2 π])))(∫_(-lim)^(+lim))(e^(-((((x-bias))^2)/(2 (σ^2)))))dx

The probability of rejection, (Subscript[P, f]), is  1 - (Subscript[P, s]).

If there are more controls, then the probability of error detection is a simple application of binomial probability

probability[bias, σ, limit, n, k] =(({{n}, {k}}))((Subscript[P, s])^k) (((1-(Subscript[P, s])))^(n-k))

where (({{n}, {k}})) is the combination of n things k at a time, n is the number of controls used, and k is the number of controls that have to be 'out' to signal a failure.

Exploring Acceptance Limits

Mathematica 6 provides the means to explore the error detection rules. The goal is to minimize false rejection — the difference between the curve and the ideal function for bias ≤ critical — while also minimizing false acceptance — the difference between the curve and the ideal function for bias > critical.

Manipulate[errorDetCurveTEERF[nConsTot,nConsOut,limit, procSD,1.0, sigmMult],
{{nConsTot, 3, "Controls"}, {3, 4, 5, 6, 9, 12, 15},ControlType -> SetterBar},
{{nConsOut, 1, "‘Out’"}, {1,2}, ControlType -> SetterBar},
{{limit,0.40, "(Limit/(Subscript[TE, a]))"}, .05, 2.0,   Appearance -> "Labeled"},
{{procSD, .20, "(SD/(Subscript[TE, a]))"}, 0.1,0.5,   Appearance -> "Labeled"},
{{sigmMult, 2.0, "(Subscript[σ, m])"}, 1.65, 3.0,  Appearance -> "Labeled"},ControlPlacement -> Left,

SaveDefinitions -> True]

Lessons learned

For any given limit, using more controls increases false rejection.

The ratio of process SD to (Subscript[TE, a])is critical.

The best rule will have a probability of error detection close to 0.9 when the bias = critical bias.

Why not let the computer compute the ideal limit?

Computing Ideal Acceptance limits

Solve the probability function P[nConsTot, nConsOut, limit, criticalBias, procSD] == errorDetectionGoal for the limit. This is easy using the Mathematica Function Solve[].

Manipulate[Quiet[QCLimitsManipulate[nTot,nOut,1.0,procSD,sm, goal]],{{nTot, 3, "Controls"}, {3,4,5,6,9, 12, 15}, ControlType-> SetterBar}, {{nOut, 1, "‘Out’"}, {1,2}}, {{procSD, 0.2, "(SD/(Subscript[TE, a]))"}, 0.099, 0.4, Appearance -> "Labeled"},
{{sm, 1.95, "(Subscript[σ, m])"}, 1.65, 3.0, Appearance -> "Labeled"},
{{goal, 0.95, "ED Goal"}, 0.5, 0.95, Appearance -> "Labeled"},ControlPlacement -> Left, SaveDefinitions -> True]

Consequences of Pass/Fail grades

Pass fail rules have too much false rejection, even when the process SD is constant. Here is how the rules work in a real situation, where the SD can vary.

Manipulate[PassFailExplorer[nTot,nOut,10.000,procSD,sm,goal],{{nTot, 3, "Controls"}, {3,4,5,6,9, 12, 15}, ControlType-> SetterBar}, {{nOut, 1, "‘Out’"}, {1,2}}, {{procSD, 2, "10\[Cross](SD/(Subscript[TE, a]))"}, 0.99, 4.0, Appearance -> "Labeled"},
{{sm, 1.95, "(Subscript[σ, m])"}, 1.65, 3.0, Appearance -> "Labeled"},
{{goal, 0.90, "ED Goal"}, 0.75, 0.95, Appearance -> "Labeled"},ControlPlacement -> Left, SaveDefinitions -> True]

Grades

Rather than just seeing if the computed Total Error ((Subscript[TE, calc])) is greater than the Allowable Total Error ((Subscript[TE, a])), the ratio is used:

r = ((Subscript[TE, calc])/(Subscript[TE, a]))

If r is much less than 1, the process performance is much better than the requirements demand. A ratio of 1.0 indicates that the process is ‘just’ acceptble, and a ratio much greater than one indicates a substantial process failure.
Using this information, a grade can be computed. Conceptually, a simple grade equation may be the following:

Grade = 100 - 30 ×((Subscript[TE, calc])/(Subscript[TE, a]))

Clearly, a grade provides much more information about the process.

Grade Computation

Manipulate[gradeVariationExplorer[pt[[1]], pt[[2]], TEa, sigMult,maxBias, nCons ], 
{{nCons, 6, "Controls"},{3, 6, 9,12, 16, 20}, ControlType -> SetterBar}, {{sigMult, 2.5, "(Subscript[σ, m])"}, 1.65, 3.0, Appearance -> "Labeled"}, {{TEa, 8.0, "(Subscript[TE, a])"}, 1, 9.5, Appearance -> "Labeled"}, {{maxBias, 0.85, "(Subscript[Bias, max])"}, 0.25, 1.0, Appearance -> "Labeled"},{{pt , {2,1.5}}, {-10, 0.05}, {10,6}, Locator},ControlPlacement -> Left, SaveDefinitions -> True]

Lessons Learned: Design

Using Mathematica, I was able to design a system so that process evaluation was clear and intuitive. Combined with an analysis of patient results, it is easy to tell the difference between a  good process and a poor one.

Using grades has proven to be a powerful means of expression, and has been used for other types of analysis: patient analysis, calibration analysis.It is far superior to pass/fail whenever the underlying model is continuous.

Good Process

Bad Process

Lessons Learned: Presentation

If a picture is worth 1000 words, an interactive graphic is worth 1000 pictures.

Communication of mathematical ideas to non-technical people is greatly enhance by using the interactive manipulations in Mathematica 6.

For most non-real-time work, programming that used to require c or c++ and specialized GUI development can be done more quickly and effectively in Mathematica.

Next Steps

All the functionality currently written in c for the production system is being redone in Mathematica 6. The resulting program will meet two needs:
1.	Provide for analysis for instruments that do not have the production system.
2.	Provide a flexible, easily adapted research tool for developing new algorithms and protocols.

Implementation

Computations and formulæ:

gauss[x_,μ_, σ_] :=Exp[-(((((x - μ)/σ))^2)/2)]/(σ(Sqrt[2 π]));

getBiasMax[TEa_, procSD_, sigmaMult_] := TEa - sigmaMult * procSD;

gaussDist[mu_, sigma_, x_] := (1/(Sqrt[2 Pi] * sigma)) E ^(-(((x - mu)/sigma)^2)/2);

probSuc1ConERF[actionLimit_, bias_, procSD_] :=(1/2) (Erf[((actionLimit-bias)/((Sqrt[2]) procSD))]+Erf[((actionLimit+bias)/((Sqrt[2]) procSD))]);

probErrDetectERF[nConsTot_, nConsOut_, actionLimit_, bias_, procSD_] :=
Module[{i, ps1},
ps1 = probSuc1ConERF[actionLimit, bias, procSD];
If[ps1<0.000001,ps1=0.000001];
If[ps1>0.999999,ps1=0.999999];
If[(nConsTot ≥ nConsOut), 1- Sum[Binomial[nConsTot, i] * ((1- ps1)^i) * (ps1 ^ (nConsTot - i)), {i, 0, nConsOut - 1}], 0]];

findAcceptanceLimitsERF[nConsTot_, nConsOut_, TEa_, procSD_, sigmaMult_,targetProbRejection_ ] :=
Module[{x, ans,maxBias},
maxBias = getBiasMax[TEa, procSD, sigmaMult];
ans = FindRoot[probErrDetectERF[nConsTot,nConsOut, x,maxBias, procSD] == targetProbRejection, {x, maxBias}, MaxIterations -> 25];
ans[[1,2]]];

falseRejectCurveTEERF[nConsTot_,nConsOut_,limit_, procSD_,TEa_]:= 
Module[{ans},
Off[NIntegrate::inum];
ans =NIntegrate[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,0,TEa/4}]/(TEa/4);
On[NIntegrate::inum];
ans
];

findBestAccLimitParametersERF[nConsTot_, nConsOut_, theTEa_, procSD_, sigmaMult_, targetProbReject_] :=
Module[{limit, falseRejection, critBias, theQCRule, ruleString},
limit = findAcceptanceLimits[nConsTot, nConsOut, theTEa, procSD, sigmaMult,targetProbReject];
falseRejection = falseRejectCurveTEERF[nConsTot,nConsOut,limit, procSD,theTEa];
critBias = getBiasMax[theTEa, procSD, sigmaMult];
theQCRule= limit/procSD;
ruleString = StringJoin["Fail if ",ToString[nConsOut], If[nConsOut == 1, " con is beyond "," cons are beyond "],ToString[theQCRule], " s, (N = ",ToString[nConsTot], ")"];
{limit, falseRejection,critBias,theQCRule, ruleString}
];

tenths[xmin_,xmax_]:=Range[Floor[xmin],Floor[xmax],0.1];
fiveTenths[xmin_,xmax_]:=Range[Floor[xmin],Floor[xmax],0.5];

errorDetCurveRuleERF[nConsTot_,nConsOut_,ruleVal_, procSD_]:=Plot[probErrDetectERF[nConsTot,nConsOut,ruleVal* procSD, x , procSD],{x,0,7* procSD},GridLines->{fiveTenths,tenths},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,ImageSize-> 500,FrameLabel->{"Bias","Prob for Rejection",StringJoin["Error Detection Graph : Fail if ",ToString[nConsOut]," con is beyond ",ToString[ruleVal]," s , N = ",ToString[nConsTot]]," "}];

errorDetCurveTEERF[nConsTot_,nConsOut_,limit_, procSD_,TEa_, sigmMult_]:=
Module[{fr,critBias, p1, p2, p3, p4},
fr = TEa/4;

critBias= getBiasMax[TEa, procSD, sigmMult];
If[fr > critBias, fr = critBias];

p1 =Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,0,fr},Filling -> Axis, FillingStyle-> Hue[0.99],GridLines->{tenths,tenths},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,ImageSize-> 433,FrameLabel->{"(Bias/(Subscript[TE, a]))","Prob for Rejection",StringJoin["Rule: Fail if ",ToString[nConsOut],If[nConsOut == 1," con is beyond ", " cons are beyond "],ToString[NumberForm[limit/procSD, {4,2}]]," s, (N = ",ToString[nConsTot],") Acceptance Limit = ±",ToString[NumberForm[limit, {4,2}]]]," "}, LabelStyle -> Directive[FontSize -> 13]];
p2 = Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,fr,TEa},GridLines->{tenths,tenths},PlotStyle -> {Blue, Thick},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,ImageSize-> 433];
p3 = Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,critBias,TEa},GridLines->{tenths,tenths},AxesOrigin->{0,0},RotateLabel->True, Filling -> 1.0, FillingStyle -> Magenta,PlotRange->All,Frame->True,ImageSize-> 433];
p4 =Graphics[{Thick, Orange, Line[{{0,0.0}, {critBias,0.0},{critBias, 1.0},{TEa, 1.0}}]}];
Show[p1, p2, p3, p4]
];

findAcceptanceLimits[nConsTot_,nConsOut_,TEa_,procSD_,sigmaMult_,targetProbRejection_]:=Module[{slope,intercept,lim1,psd1,lim2,psd2,ans},psd1=0.01*TEa;
lim1=findAcceptanceLimitsERF[nConsTot,nConsOut,TEa,psd1,sigmaMult,targetProbRejection];
psd2=0.15*TEa;
lim2=findAcceptanceLimitsERF[nConsTot,nConsOut,TEa,psd2,sigmaMult,targetProbRejection];
slope=(lim2-lim1)/(psd2-psd1);
intercept=lim1-slope*psd1;
ans=procSD*slope+intercept;
ans=If[(ans>0),ans,0];
ans
]

errorDetCurve[nConsTot_,nConsOut_,limit_, procSD_,TEa_, sigmMult_]:=
Module[{fr,critBias, p1, p2, p3, p4, p5},
fr = TEa/4;
critBias= getBiasMax[TEa, procSD, sigmMult];
If[fr > critBias, fr = critBias];

p1 =Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,0,fr},Filling -> Axis, FillingStyle-> Hue[0.99],GridLines->{tenths,tenths},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,LabelStyle  ->Directive[FontSize->13],ImageSize-> 433,FrameLabel->{"(Bias/(Subscript[TE, a]))","Prob for Rejection",StringJoin["Rule: Fail if ",ToString[nConsOut],If[nConsOut == 1," con is beyond ", " cons are beyond "],ToString[NumberForm[limit/procSD, {4,2}]]," s, (N = ",ToString[nConsTot],") Acceptance Limit = ±",ToString[NumberForm[limit, {4,2}]]]," "}];
p2 = Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,fr,TEa},GridLines->{fiveTenths,tenths},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,ImageSize-> 433];
p3 = Plot[probErrDetectERF[nConsTot,nConsOut,limit, x , procSD],{x,critBias,TEa},GridLines->{fiveTenths,tenths},AxesOrigin->{0,0},RotateLabel->True,PlotRange->All,Frame->True,ImageSize-> 433, Filling -> Top, FillingStyle -> Magenta];
p4 =Graphics[{Thick, Orange, Line[{{0,0.0}, {critBias,0.0},{critBias, 1.0},{TEa, 1.0}}]}];
p5 =Graphics[{Thick, Green, Line[{{limit,0.5},{limit, 1.0}}]}];
Show[ p1, p2,p3,  p4, p5]
];

QCLimitsManipulate[nConsTot_,nConsOut_,TotErr_,procSD_,sigmaMult_, targetProbReject_] :=
Module[{parameters, limit, critBias, falseAcc},
critBias= getBiasMax[TotErr, procSD, sigmaMult];
parameters =findBestAccLimitParametersERF[nConsTot,nConsOut,TotErr,procSD,sigmaMult, targetProbReject];
limit = parameters[[1]];
falseAcc = falseAccCurveTEERF[nConsTot,nConsOut,critBias, limit,procSD,TotErr];
Column[{Row[{"(Subscript[TE, a]) set to 1.0."}],
Row[{"Green   : acceptance limit = ±", NumberForm[ parameters[[1]], {4,2}]}],Row[{"Orange  : ideal ED, crit. bias = ",NumberForm[ critBias, {4,2}]}],Row[{"Red     : ", NumberForm[100.0  parameters[[2]], {4,2}], "% false rejection."}],
Row[{"Magenta : ", NumberForm[100 *  falseAcc, {5,3}], "% false acceptance."}],
errorDetCurve[nConsTot,nConsOut,limit, procSD,TotErr, sigmaMult]}]

];

PassFailExplorer[nConsTot_,nConsOut_,TotErr_,procSD_,sigmaMult_, targetProbReject_] :=
Module[{parameters, limit, critBias, falseAcc, p1, p2},
critBias= getBiasMax[TotErr, procSD, sigmaMult];
parameters =findBestAccLimitParametersERF[nConsTot,nConsOut,TotErr,procSD,sigmaMult, targetProbReject];
limit = parameters[[1]];
p1 = ContourPlot[probErrDetectERF[nConsTot, nConsOut, limit, bias, SD], {bias, 0, TotErr}, {SD, 0, TotErr}, Contours -> {0.25, 0.5, 0.75, targetProbReject},Mesh -> None, ContourShading -> {None, None, None, None,Opacity[.60, LightRed]}, ContourStyle->{Directive[ Thick, Blue],Directive[ Thick, Orange], Directive[ Thick, Magenta], Directive[ Thick,Red]}, ImageSize ->  433];

p2 = Plot[-x/sigmaMult + TotErr/sigmaMult, {x, 0, TotErr},Filling -> Axis, FillingStyle -> RGBColor[0,1,0, .5], PlotRange -> {{0, TotErr}, {0, TotErr}}, AspectRatio -> Automatic, Frame -> True, GridLines -> Automatic, FrameLabel -> {"| bias |","true SD"," "," "}, ImageSize ->  433]; 
Show[p2, p1]

]

conGrade2[bias_, sd_, sm_, tea_, maxBias_] :=
Module[{x1, y1, m, biasCutoff, r, grade, absBias},
biasCutoff = maxBias * tea;
absBias = Abs[bias];
m = 1.0;
If[ (absBias == 0.0), {x1 = 0.0; y1 = tea/sm},
{m =  sd/absBias;
x1 = tea/(sm * m +1);
y1 = tea * m/(sm * m + 1)}];
If[x1 > biasCutoff, {x1 = biasCutoff; y1 = m * biasCutoff}];
r = (absBias + sd)/(Abs[x1] + Abs[y1]);
grade = 104.0 - 34.0 r;
If[grade > 100.0, grade = 100.0];
If[grade < 0.0, grade = 0.0];
grade
]

testFuncSimple[ntrials_, testMean_, testSD_,sm_, tea_, maxBias_] :=
 Module[{myList, compMean, compSD, grade},
myList = RandomReal[NormalDistribution[testMean,testSD], ntrials];
compMean = Mean[myList];
compSD = StandardDeviation[myList];
grade = conGrade2[compMean, compSD, sm,tea, maxBias];
{grade, {compMean,compSD}}
]

Needs["Histograms`"]

gradeDistGraph[μ_, σ_, theData_List] := Module[{curve, hist, scale, maxBin},
maxBin = Max[BinCounts[theData,{0, 100, 1}]];
scale = maxBin/gauss[μ, μ, σ];
curve = Plot[ scale gauss[x, μ, σ], {x, 0, 100}, PlotStyle -> Thick, PlotRange -> {{0, 100}, All}];
hist = Histogram[theData, HistogramRange -> {0, 100}, HistogramCategories -> 100, GridLines -> Automatic];
Show[hist,curve, PlotRange -> All, ImageSize -> {433, 200}, AspectRatio -> 3/10, Frame -> True, FrameLabel -> {{"frequency", " "}, {"grade", "Distribution of computed grades"}}]]

gradePlotBkgd[sm_, tea_, maxBias_] := Module[{p1, myPoly},
myPoly =Polygon[{{-tea * maxBias, 0},{-tea * maxBias, (tea *(1 - maxBias))/sm}, {0, tea/sm}, {tea * maxBias, (tea *(1- maxBias))/sm},{tea * maxBias, 0}, {-tea * maxBias, 0}}]; Graphics[{EdgeForm[Directive[Thickness[0.01], Red]],Directive[ Green, Opacity[ 0.7]], myPoly}, PlotRange -> {{-10.1, 10.1}, {0, 6}}, GridLines -> Automatic, Frame ->True, LabelStyle ->Directive[FontSize -> 14],FrameLabel -> {{"S.D.", ""},{"bias", ""}}, ImageSize ->{433, 300}, AspectRatio -> Full]
]

fancyGauss[μ_, σ_, tea_, range_, vertRange_] := Module[ {p1, p2, p3, outside,close},
outside = 1.0 - NIntegrate[gauss[x, μ, σ], {x,-tea, tea}, AccuracyGoal -> 7];
close =  NIntegrate[gauss[x, μ, σ], {x,-tea/2, tea/2}, AccuracyGoal -> 7];
If[outside < .0000001, outside = 0.0];
If[close < .0000001, close = 0.0];
p1 = Plot[gauss[x, μ, σ], {x, -range, -tea}, Filling -> Axis, FillingStyle -> Red,GridLines  -> Automatic, PlotRange -> {{-range, range}, {0, vertRange}}];p2 = Plot[gauss[x, μ, σ], {x, -tea, tea}, Filling -> Axis, FillingStyle -> Green, PlotRange -> {{-range, range}, {0, vertRange}}];
p3 =Plot[gauss[x, μ, σ], {x, tea, range}, Filling -> Axis, FillingStyle -> Red, PlotRange -> {{-range, range}, {0, vertRange}}];
Show[p1,p2,p3, ImageSize -> {500, 250}, AspectRatio -> 3/10,Frame ->True, LabelStyle -> Directive[FontSize -> 14],FrameLabel -> {{"frequency", ""}, {"bias", "raw defect rate: "<>ToString[NumberForm[100 *outside, {6, 6}]]<>"%"<>"   ‘Close’ to True Val : "<> ToString[NumberForm[100 * close, {6, 6}]] <> "%"}}]]

simpleGradeExplorer[nomBias_, nomSD_, tea_, sm_,maxBias_ ] :=Module[{nomGrade, gradeMean, gradeSD, gradesList, pointsList},

Column[{
Row[{"Grade = ",NumberForm[conGrade2[nomBias, nomSD, sm, tea, maxBias], {3, 1}],"  nominal Bias = ", NumberForm[nomBias, {4,2}], "   nominal SD = ", NumberForm[nomSD, {4,2}]}],
gradePlotBkgd[sm, tea, maxBias],
fancyGauss[nomBias, nomSD, tea,10, .8]
}, 0 ] 
]

gradeTestForPlotSimple[ncontrols_, testMean_, testSD_,sm_, tea_, maxBias_, nTrials_] :=
Module[{trueGrade, rawList,trialsList,pointsList, gradeMean, gradeSD},
 trueGrade = conGrade2[testMean, testSD, sm, tea, maxBias];
rawList = Table[testFuncSimple[ncontrols, testMean, testSD, sm, tea, maxBias], {nTrials}];
trialsList = rawList[[All,1]];
pointsList = rawList[[All, 2]];
gradeMean = Mean[trialsList];
gradeSD = StandardDeviation[trialsList];
{trueGrade, gradeMean, gradeSD, trialsList, pointsList}]

gradeVariationExplorer[nomBias_, nomSD_, tea_, sm_,maxBias_, nCons_] :=Module[{nomGrade, gradeMean, gradeSD, gradesList, pointsList, p1, p2},
{nomGrade, gradeMean, gradeSD,gradesList, pointsList} = gradeTestForPlotSimple[nCons, nomBias, nomSD,sm, tea, maxBias,1000];
p1 = gradePlotBkgd[sm, tea, maxBias];
p2 = ListPlot[pointsList, PlotRange -> {{-10, 10}, {0,6}, All}];

Column[{
Row[{"Grade = ", NumberForm[conGrade2[nomBias, nomSD, sm, tea, maxBias], {3, 1}],"  Bias: ", NumberForm[nomBias, {4,2}], "  SD: ", NumberForm[nomSD, {4,2}]}],Row[{"CI based on 1000 trials."}],Row[{"Grade 95% CI = ±", NumberForm[2 gradeSD, {4,2}], " (", NumberForm[nomGrade - 2 gradeSD, {4,1}], " - ", NumberForm[nomGrade + 2 gradeSD, {4,1}], ")"}],
Show[p1, p2],
gradeDistGraph[gradeMean, gradeSD, gradesList]
}] 
]