Model Checking Diagnostics
• Assumptions
1 Model is correct
Model Adequacy 2 Independent observations
3 Errors normally distributed
4 Constant variance
Design of Experiments - Montgomery
Sections 3-4 and 14-1.1
yij = (y .. + (y i. − y .. )) + (yij − y i.)
yij = ŷij + ²̂ij
observed = predicted + residual
• Diagnostics use predicted and residual values
5 5-1
Diagnostics Normality Assumption
• Histogram
• Normality
– Is histogram of residuals bell-shaped?
– Histogram of residuals
– Normal probability plot / QQ plot • Normal Probability Plot
– Shapiro-Wilks/Kolmogorov-Smirnov Test – Ordered residuals vs cumulative probability points
• Variance – On normal probability paper
– Plot ²̂ij vs ŷij (residual plot) ∗ Pt = (t − .5)/N for t = 1, 2, . . . , N
– Bartlett’s or Levene’s Test
∗ Plot Pt vs ²̂(t)
• Independence – Is it approximately linear?
– Plot ²̂ij vs time/space
– Plot ²̂ij vs variable of interest • QQ plot
– Similar plot but compute associated quantile values
• Outliers
– Is it influential? With and without analysis – Find zt such that P(Z ≤ zt ) = Pt
– Formal tests (e.g. standardized residuals)
– On regular graph paper (zt vs ²̂(t) )
– Investigate why result may occur, don’t try to eliminate
5-2 5-3
Constant Variance
• Often experiments with non-constant variance
• Size of residual associated with predicted value Independence
• Residual plot
• Plot of the residuals over time
– Plot ²̂ij vs ŷij
– Is there a drift or pattern as trials proceed?
– Is the range constant for different levels of ŷij
• Bartlett’s Test (page 81-82) • Plot residuals versus relevant variables
– More formal test – Often variables omitted from analysis
– Compares pooled with a sample variances
– Experimental conditions (e.g., temp)
– Often unnecessary
– Very sensitive to Normality assumption – May result in inclusion of factor in next exp
• Levene’s Test
– Compute |yij − mi | where mi is median
– Compare samples using F-test
5-4 5-5
options ls=80 ps=65; Diagnostics Example
title1 ’Diagnostics Example’; The GLM Procedure
data one; Dependent Variable: strength
infile ’c:\saswork\data\[Link]’; Sum of
input percent strength time; Source DF Squares Mean Square F Value Pr > F
Model 4 475.7600000 118.9400000 14.76 <.0001
proc glm data=one; Error 20 161.2000000 8.0600000
class percent; Corrected Total 24 636.9600000
model strength=percent;
means percent / hovtest=bartlett hovtest=levene; R-Square Coeff Var Root MSE strength Mean
output out=diag p=pred r=res; 0.746923 18.87642 2.839014 15.04000
proc sort; by pred; Source DF Type I SS Mean Square F Value Pr > F
symbol1 v=circle i=sm50; title1 ’Residual Plot’; percent 4 475.7600000 118.9400000 14.76 <.0001
proc gplot; plot res*pred/frame; run;
Source DF Type III SS Mean Square F Value Pr > F
proc univariate data=diag noprint; percent 4 475.7600000 118.9400000 14.76 <.0001
var res; qqplot res / normal (L=1 mu=est sigma=est);
histogram res / normal; run; Levene’s Test for Homogeneity of strength Variance
run; ANOVA of Squared Deviations from Group Means
proc sort; by time; Sum of Mean
symbol1 v=circle i=sm75; Source DF Squares Square F Value Pr > F
title1 ’Plot of residuals vs time’; percent 4 91.6224 22.9056 0.45 0.7704
proc gplot; plot res*time / vref=0 vaxis=-6 to 6 by 1; Error 20 1015.4 50.7720
run;
symbol1 v=circle i=sm50; Bartlett’s Test for Homogeneity of strength Variance
title1 ’Plot of residuals vs time’;
proc gplot; plot res*time / vref=0 vaxis=-6 to 6 by 1; Source DF Chi-Square Pr > ChiSq
run; percent 4 0.9331 0.9198
5-6 5-7
Constant Variance
• Often experiments with non-constant variance
• Does not affect F-test dramatically (balanced)
• Why concern?
– Comparison of treatments depends on MSE
– Incorrect intervals and comparison results
• Variance-Stabilizing Transformations
– Common transformations
√ √ √
x, log(x), 1/x, arcsin( x), and 1/ x
– Box-Cox transformations
uses maximum likelihood procedures
can approximate using relationship σi = θµβi
transformation is X 1−β
• Dist often more “normal” after transformation
5-8 5-9
Transformations
• Consider response X with E(X)=µx and Var(X)=σx2 Box-Cox Transformation
• Define Y = f (X); What is the mean and var of Y ?
• Perform analysis of variance on
Delta Method
y λ −1
Consider f (X) where f 0 (µx ) 6= 0 λ 6= 0
λ λẏ λ−1
f (X) ≈ f (µx ) + (X − µx )f 0 (µ x) y =
ẏlogy λ = 0
E(Y )=E(f (X))≈ E(f (µx )) + E((X − µx )f 0 (µx ))= f (µx )
Var(Y ) ≈ [f 0 (µx )]2 Var(X) = [f 0 (µx )]2 σx2 ẏ is the geometric mean of the observations
• Suppose σx2 depends on µx → σx2 = g(µx) 1/N
Y
a Y
ni
• Want to find Y = f (X) such that Var(Y )≈ c ẏ = yij
i=1 j=1
• Have shown Var(f (X))≈ [f 0 (µx)]2 σx2
• y to the power of λ rescaled for direct comparison
• Want to choose f such that [f 0 (µx )]2g(µx) ≈ c
• Find λ which minimizes SSE
Examples
R √
g(µ) = µ (Poisson) f (X) =
R
√1 dµ → f (X) = X • [Link] does these calculations for you
µ √
g(µ) = µ(1 − µ) (Binomial) f (X) = √ 1 dµ → f (X) = asin( X)
R µ(1−µ)
g(µ) = µ2β (Box-Cox) f (X) = R µ−β dµ → f (X) = X 1−β
g(µ) = µ2 (Box-Cox) f (X) = 1
µ
dµ → f (X) = log X
5-10 5-11
[Link]
options nocenter ps=65 ls=80;
title1 ’Increasing Variance Example’;
data one;
infile ’c:\saswork\data\[Link]’; input trt resp;
proc glm data=one; class trt;
model resp=trt; output out=diag p=pred r=res;
title1 ’Residual Plot’; symbol1 v=circle i=none;
proc gplot data=diag; plot res*pred /frame;
proc univariate data=one noprint;
var resp; by trt; output out=two mean=mu std=sigma;
data three;
set two; logmu = log(mu); logsig = log(sigma);
proc reg; model logsig = logmu;
title1 ’Mean vs Std Dev’; symbol1 v=circle i=rl;
proc gplot; plot logsig*logmu / regeqn;
proc transreg data=one;
model boxcox(resp / lambda=-2 to 2 by .2) = class(trt);
run;
5-12 5-13
Transformation Information for BoxCox(resp)
Lambda R-Square Log Like OBS L SS OBS L SS
1 -2.00 2150.06 10 0.25 112.37
-2.0 0.10 -108.906 2 -1.75 1134.83 11 0.50 154.23
-1.8 0.10 -93.460 3 -1.50 628.94 12 0.75 253.63
-1.6 0.11 -78.920 4 -1.25 369.35 13 1.00 490.36
-1.4 0.12 -65.416 5 -1.00 232.32 14 1.25 1081.29
-1.2 0.13 -53.103 6 -0.75 158.56 15 1.50 2636.06
7 -0.50 119.28 16 1.75 6924.95
-1.0 0.14 -42.153 8 -0.25 100.86 17 2.00 19233.39
-0.8 0.16 -32.767 9 0.00 98.09
-0.6 0.17 -25.184
-0.4 0.19 -19.683
-0.2 0.20 -16.593 *
0.0 + 0.21 -16.284 <
0.2 0.22 -19.124
0.4 0.22 -25.398
0.6 0.22 -35.200
0.8 0.20 -48.378
1.0 0.18 -64.563
1.2 0.17 -83.272
1.4 0.15 -104.019
1.6 0.13 -126.377
1.8 0.12 -150.006
2.0 0.10 -174.641
< - Best Lambda
* - Confidence Interval
+ - Convenient Lambda
5-14 5-15
Unusual Observations Near-Zero/Truncated Values
• Can use residual plot to detect unusual observations • Can have heavily skewed observations (near zero)
– Concentration of rare contaminant
• Unusual = large |²̂ij |
– Number of defects in assembly line
• Sometimes typographical error – Number of birds at a given site
• Otherwise worth investigation • Measurements may be truncated
– Helpful to have detailed lab notes – Concentration of rare contaminant (not detectable)
• Don’t immediately look for excuses – Lifetime of component (does not fail)
• Transformations sometimes successful
√
• Formal test to see if unusual log(x+.001), x + .001
• But does not answer if in error • Non-parametric analysis (Kruskal-Wallis)
• More advanced
Assume non-Normal errors (GLIM models)
• Does observation make a difference if excluded?
EM algorithm (iterative MLE approach)
• Perform analysis with and without observations(s) Bayesian methods
5-16 5-17
options nocenter ps=65 ls=80; N P A R 1 W A Y P R O C E D U R E
data new; Wilcoxon Scores (Rank Sums) for Variable NITROGEN
input strain nitrogen @@; Sum of Expected Std Dev Mean
cards; STRAIN N Scores Under H0 Under H0 Score
1 19.4 1 32.6 1 27.0 1 32.1 1 33.0
2 17.7 2 24.8 2 27.9 2 25.2 2 24.3 1 5 130.0 77.5000000 17.9618849 26.0000000
3 17.0 3 19.4 3 9.1 3 11.9 3 15.8 2 5 111.0 77.5000000 17.9618849 22.2000000
4 20.7 4 21.0 4 20.5 4 18.8 4 18.6 3 5 40.0 77.5000000 17.9618849 8.0000000
5 14.3 5 14.4 5 11.8 5 11.6 5 14.2 4 5 88.0 77.5000000 17.9618849 17.6000000
6 17.3 6 19.4 6 19.1 6 16.9 6 20.8 5 5 23.0 77.5000000 17.9618849 4.6000000
; 6 5 73.0 77.5000000 17.9618849 14.6000000
Average Scores Were Used for Ties
proc npar1way;
class strain; Kruskal-Wallis Test (Chi-Square Approximation)
var nitrogen; CHISQ = 21.659 DF = 5 Prob > CHISQ = 0.0006
run; Median Scores (Number of Points Above Median) for Variable NITROGEN
___________________________________________________
Sum of Expected Std Dev Mean
N P A R 1 W A Y P R O C E D U R E STRAIN N Scores Under H0 Under H0 Score
Analysis of Variance for Variable NITROGEN 1 5 5.0 2.50000000 1.03806850 1.00000000
2 5 4.0 2.50000000 1.03806850 0.80000000
STRAIN N Mean Among MS Within MS 3 5 1.0 2.50000000 1.03806850 0.20000000
169.409333 11.7886667 4 5 3.0 2.50000000 1.03806850 0.60000000
1 5 28.8200000 5 5 0.0 2.50000000 1.03806850 0.00000000
2 5 23.9800000 F Value Prob > F 6 5 2.0 2.50000000 1.03806850 0.40000000
3 5 14.6400000 14.371 0.0001 Average Scores Were Used for Ties
4 5 19.9200000
5 5 13.2600000 Median 1-Way Analysis (Chi-Square Approximation)
6 5 18.7000000 CHISQ = 13.533 DF = 5 Prob > CHISQ = 0.0189
5-18 5-19