Stat 512 Class 1 - Purdue University

Download Report

Transcript Stat 512 Class 1 - Purdue University

Topic 31: Two-way
Random Effects Models
Outline
• Two-way random effects design
– Model
– F tests
Data for two-way design
• Y is the response variable
• Factor A with levels i = 1 to a
• Factor B with levels j = 1 to b
• Yijk is the kth observation in cell (i, j)
k = 1 to nij
• Have balanced designs with n = nij
KNNL Example
• KNNL Problem 25.15, p 1080
• Y is fuel efficiency in miles per gallon
• Factor A represents four different
drivers, a=4 levels
• Factor B represents five different cars
of the same model , b=5
• Each driver drove each car twice over
the same 40-mile test course
Read and check the data
data a1;
infile 'c:\...\CH25PR15.TXT';
input mpg driver car;
proc print data=a1;
run;
The data
Obs
1
2
3
4
5
6
7
8
9
10
mpg
25.3
25.2
28.9
30.0
24.8
25.1
28.4
27.9
27.1
26.6
driver
1
1
1
1
1
1
1
1
1
1
car
1
1
2
2
3
3
4
4
5
5
Prepare the data for a plot
data a1; set a1;
if (driver eq 1)*(car eq 1)
then dc='01_1A';
if (driver eq 1)*(car eq 2)
then dc='02_1B';
⋮
if (driver eq 4)*(car eq 5)
then dc='20_4E';
Plot the data
title1 'Plot of the data';
symbol1 v=circle i=none c=black;
proc gplot data=a1;
plot mpg*dc/frame;
run;
Find the means
proc means data=a1;
output out=a2 mean=avmpg;
var mpg;
by driver car;
Plot the means
title1 'Plot of the means';
symbol1 v='A' i=join c=black;
symbol2 v='B' i=join c=black;
symbol3 v='C' i=join c=black;
symbol4 v='D' i=join c=black;
symbol5 v='E' i=join c=black;
proc gplot data=a2;
plot avmpg*driver=car/frame;
run;
Random effects model
• Yijk = μij + εijk
– the μij are N(μ, σμ2)
– the εijk are iid N(0, σ2)
– μij and εijk are independent
• Dependence among the Yijk can be
most easily described by specifying
the covariance matrix of the vector
(Yijk)’
Random factor effects
model
• Yijk = μ + i + j + ()ij + εijk
• i ~ N(0, σ2)
• j ~ N(0, σ2)
• ()ij ~ N(0, σ2)
• εij ~ N(0, σ2)
• σY2 = σ2 + σ2 + σ2 + σ2
Parameters
• There are five parameters in this
model
–μ
– σ 2
– σ 2
– σ2
– σ2
• The cell means are random variables,
not parameters
ANOVA table
• The terms and layout of the ANOVA
table are the same as what we used
for the fixed effects two-way model
• The expected mean squares (EMS)
are different
EMS and parameter
estimates
•
•
•
•
•
E(MSA) = σ2 + bnσ2 + nσ2
E(MSB) = σ2 + anσ2 + nσ2
E(MSAB) = σ2 + nσ2
E(MSE) = σ2
Estimates of the variance components
can be obtained from these equations,
replacing E(MS) with table value, or other
methods such as ML
Hypotheses
• H0A: σ2 = 0; H1A: σ2 ≠ 0
• H0B: σ2 = 0; H1B : σ2 ≠ 0
• H0AB : σ2 = 0; H1AB : σ2 ≠ 0
Hypothesis H0A
•
•
•
•
•
H0A: σ2 = 0; H1A: σ2 ≠ 0
E(MSA) = σ2 + bnσ2 + nσ 2
E(MSE) = σ2
E(MSAB) = σ2 + nσ 2
So H0A is tested by F = MSA/MSAB
with df a-1 and (a-1)(b-1)
Hypothesis H0B
•
•
•
•
H0B: σ2 = 0; H1B: σ2 ≠ 0
E(MSB) = σ2 + anσ2 + nσ 2
E(MSAB) = σ2 + nσ 2
So H0B is tested by F = MSB/MSAB
with df b-1 and (a-1)(b-1)
Hypothesis H0AB
•
•
•
•
H0AB: σ 2 = 0; H1AB: σ2 ≠ 0
E(MSAB) = σ2 + nσ2
E(MSE) = σ2
So H0AB is tested by F = MSAB/MSE
with df (a-1)(b-1) and ab(n-1)
Run proc glm
proc glm data=a1;
class driver car;
model mpg=driver car
driver*car;
random driver car
driver*car/
test;
run;
Model and error output
Source
Model
Sum of
DF
Squares
19 377.444750
Error
20 3.5150000
Corrected Total 39 380.959750
Mean
Square F Value Pr > F
19.8655132 113.03 <.0001
0.1757500
Factor effects output
All using MSE in denominator
Source
driver
car
DF
Type I SS Mean Square F Value Pr > F
3 280.2847500
93.4282500 531.60 <.0001
4 94.7135000
driver*car
12
2.4465000
23.6783750
134.73 <.0001
0.2038750
1.16 0.3715
Only the interaction test
is valid here
Random statement output
Source
driver
Type III Expected Mean Square
Var(Error) + 2 Var(driver*car) + 10 Var(driver)
car
Var(Error) + 2 Var(driver*car) + 8 Var(car)
driver*car
Var(Error) + 2 Var(driver*car)
Random/test output
Source
driver
car
Error
DF Type III SS Mean Square F Value Pr > F
3 280.284750
93.428250 458.26 <.0001
4
94.713500
23.678375
12
2.446500
0.203875
Error: MS(driver*car)
116.14 <.0001
Random/test output
Source
driver*car
Error:
MS(Error)
DF Type III SS
12 2.446500
20
3.515000
Mean
Square F Value Pr > F
0.203875
1.16 0.3715
0.175750
Proc varcomp
proc varcomp data=a1;
class driver car;
model mpg=driver car
driver*car;
run;
Output
MIVQUE(0) Estimates
Variance Component
Var(driver)
mpg
9.32244
Var(car)
2.93431
Var(driver*car)
0.01406
Var(Error)
0.17575
MIXED Output
Fit Statistics
-2 Res Log Likelihood
86.8
AIC (smaller is better)
94.8
AICC (smaller is better) 96.0
BIC (smaller is better)
93.2
Covariance Parameter
Estimates
Cov Parm
Estimate
car
2.9343
driver
9.3224
car*driver
0.01406
Residual
0.1757
DEFAULT: F-tests only supplied
for fixed effects
Estimated V Correlation Matrix for Subject 1
Row Col1 Col2 Col3 Col4 Col5 Col6 Col7 Col8 Col9 Col10 Col11 Col12 Col13 Col14 Col15
1 1.0000 0.9859 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.2358 0.2358
2 0.9859 1.0000 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490 0.2358 0.2358
3 0.7490 0.7490 1.0000 0.9859 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490
0.2358 0.2358
4 0.7490 0.7490 0.9859 1.0000 0.7490 0.7490 0.7490 0.7490 0.7490 0.7490
0.2358 0.2358
5 0.7490 0.7490 0.7490 0.7490 1.0000 0.9859 0.7490 0.7490 0.7490 0.7490
0.2358
6 0.7490 0.7490 0.7490 0.7490 0.9859 1.0000 0.7490 0.7490 0.7490 0.7490
0.2358
Covariances/Correlations
Var (Yijk )             12.447
2
2
2
2
Cov (Yijk , Yijk  )           12.271
2
2
2
Corr (Yijk , Yijk  )  12.271 / 12.447  0.9859
Cov (Yijk , Yijk  )     9.322
2
Corr (Yijk , Yijk  )  9.322 / 12.447  0.7489
Using Proc Mixed
proc mixed data=a1 covtest;
class car driver;
model mpg=;
random car driver
car*driver/vcorr;
run;
MIXED Output
Covariance Parameter Estimates
Standard
Z
Cov Parm Estimate
Error Value Pr > Z
car
2.9343
2.0929 1.40 0.0805
driver
car*driver
Residual
9.3224
7.6284
1.22 0.1108
0.01406
0.05004
0.28 0.3893
0.1757
0.05558
3.16 0.0008
Last slide
• Finish reading KNNL Chapter 25
• We used program topic31.sas to
generate the output for today