The BMI Health study with SAS proc calis
/********************** bmi1.sas **************************/
title 'BMI and Health: Read data and analyze ignoring measurement error';
data health;
infile '/folders/myfolders/431s15/bmihealth.data.txt';
input age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
/* fat1 and fat2 are percent body fat */
age = (age1+age2)/2; bmi = (bmi1+bmi2)/2; fat = (fat1+fat2)/2;
cholest = (cholest1+cholest2)/2 ; diastol = (diastol1+diastol2)/2;
proc means;
var age1 -- diastol;
ods graphics off; /* Suppress regression plots */
proc reg;
title2 'Regression on average measurements';
model cholest diastol = age bmi fat;
BMI: mtest bmi=0; /* Multivariate test */
/********************** bmi2.sas **************************/
title 'BMI and Health: Use the Double Measurement Design';
data health;
infile '/folders/myfolders/431s15/bmihealth.data.txt';
input age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
/* fat1 and fat2 are percent body fat */
age = (age1+age2)/2; bmi = (bmi1+bmi2)/2; fat = (fat1+fat2)/2;
cholest = (cholest1+cholest2)/2 ; diastol = (diastol1+diastol2)/2;
proc calis pshort nostand vardef=n pcorr;
/* pshort and nostand suppress some output.
vardef=n uses n in the denominator of Sigmahat
rather than the default n-1, and makes Chisquare
for model fit a traditional likelihood ratio test.
pcorr prints correlation (actually, covariance)
matrices of observed variables -- both sample and
reproduced. */
title2 'Full Model';
var /* Name the observed variables */
age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
/* Now give simultaneous equations, separated by commas. Latent
variables begin with F for factor. Error terms begin with
E for error or D for disturbance. SAS is not case sensitive.
You must name all the parameters. Optional starting values in
parentheses may be given after the parameters. */
lineqs
Fcholest = beta11 Fage + beta12 Fbmi + beta13 Ffat + epsilon1,
Fdiastol = beta21 Fage + beta22 Fbmi + beta23 Ffat + epsilon2,
age1 = Fage + e11, bmi1 = Fbmi + e12,
fat1 = Ffat + e13,
cholest1 = Fcholest + e14, diastol1 = Fdiastol + e15,
age2 = Fage + e21,bmi2 = Fbmi + e22,
fat2 = Ffat + e23,
cholest2 = Fcholest + e24, diastol2 = Fdiastol + e25;
variance /* Variances of exogenous vars */
Fage = phi11, Fbmi = phi22, Ffat = phi33,
epsilon1 = psi11, epsilon2 = psi22,
e11 = omega111, e12 = omega122, e13 = omega133,
e14 = omega144, e15 = omega155,
e21 = omega211, e22 = omega222, e23 = omega233,
e24 = omega244, e25 = omega255;
cov /* Covariances */
Fage Fbmi = phi12, Fage Ffat = phi13, Fbmi Ffat = phi23,
epsilon1 epsilon2 = psi12,
e11 e12 = omega112, e11 e13 = omega113, e11 e14 = omega114,
e11 e15 = omega115,
e12 e13 = omega123, e12 e14 = omega124, e12 e15 = omega125,
e13 e14 = omega134, e13 e15 = omega135,
e14 e15 = omega145,
e21 e22 = omega212, e21 e23 = omega213, e21 e24 = omega214,
e21 e25 = omega215,
e22 e23 = omega223, e22 e24 = omega224, e22 e25 = omega225,
e23 e24 = omega234, e23 e25 = omega235,
e24 e25 = omega245;
bounds /* Variances are positive */
0.0 < phi11 phi22 phi33 psi11 psi22
omega111 omega122 omega133 omega144 omega155
omega211 omega222 omega233 omega244 omega255;
/* Now fit a reduced model to test H0: beta12 = beta22 = 0,
meaning BMI is unrelated to either cholesterol or blood pressure
if we allow for age and percent body fat. Copy the code; only
the last line is different. */
proc calis pshort nostand vardef=n pcorr;
title2 'Reduced Model with beta12=beta22=0';
var /* Name the observed variables */
age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
/* Now give simultaneous equations, separated by commas. Latent
variables begin with F for factor. Error terms begin with
E for error or D for disturbance. SAS is not case sensitive.
You must name all the parameters. Optional starting values in
parentheses may be given after the parameters. */
lineqs
Fcholest = beta11 Fage + beta12 Fbmi + beta13 Ffat + epsilon1,
Fdiastol = beta21 Fage + beta22 Fbmi + beta23 Ffat + epsilon2,
age1 = Fage + e11, bmi1 = Fbmi + e12,
fat1 = Ffat + e13,
cholest1 = Fcholest + e14, diastol1 = Fdiastol + e15,
age2 = Fage + e21,bmi2 = Fbmi + e22,
fat2 = Ffat + e23,
cholest2 = Fcholest + e24, diastol2 = Fdiastol + e25;
variance /* Variances of exogenous vars */
Fage = phi11, Fbmi = phi22, Ffat = phi33,
epsilon1 = psi11, epsilon2 = psi22,
e11 = omega111, e12 = omega122, e13 = omega133,
e14 = omega144, e15 = omega155,
e21 = omega211, e22 = omega222, e23 = omega233,
e24 = omega244, e25 = omega255;
cov /* Covariances */
Fage Fbmi = phi12, Fage Ffat = phi13, Fbmi Ffat = phi23,
epsilon1 epsilon2 = psi12,
e11 e12 = omega112, e11 e13 = omega113, e11 e14 = omega114,
e11 e15 = omega115,
e12 e13 = omega123, e12 e14 = omega124, e12 e15 = omega125,
e13 e14 = omega134, e13 e15 = omega135,
e14 e15 = omega145,
e21 e22 = omega212, e21 e23 = omega213, e21 e24 = omega214,
e21 e25 = omega215,
e22 e23 = omega223, e22 e24 = omega224, e22 e25 = omega225,
e23 e24 = omega234, e23 e25 = omega235,
e24 e25 = omega245;
bounds /* Variances are positive */
0.0 < phi11 phi22 phi33 psi11 psi22
omega111 omega122 omega133 omega144 omega155
omega211 omega222 omega233 omega244 omega255;
lincon beta12=0, beta22=0;
proc iml;
title2 'Calculate Likelihood ratio test of H0: beta12=beta22=0';
G2 = 500 * (0.0122914464 - 0.0093074908);
/* Difference between final objective function values */
pval = 1 - probchi(G2,2);
print G2 pval;
print "Or better, difference between 'Chi-Square' values";
diff = 6.1457-4.6537;
print "6.1457-4.6537 = " diff;
/********************** bmi3.sas **************************/
title 'BMI and Health: Like bmi2.sas, but try to make it shorter';
title2 'Also Try Wald Test: Compare G2 = 1.492, df=2, p = 0.474';
data health;
infile '/folders/myfolders/431s15/bmihealth.data.txt';
input age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
/* fat1 and fat2 are percent body fat */
age = (age1+age2)/2; bmi = (bmi1+bmi2)/2; fat = (fat1+fat2)/2;
cholest = (cholest1+cholest2)/2 ; diastol = (diastol1+diastol2)/2;
proc calis pshort nostand vardef=n pcorr;
title3 'Full Model Only';
var /* Name the observed variables */
age1 bmi1 fat1 cholest1 diastol1
age2 bmi2 fat2 cholest2 diastol2;
lineqs
Fcholest = beta11 Fage + beta12 Fbmi + beta13 Ffat + epsilon1,
Fdiastol = beta21 Fage + beta22 Fbmi + beta23 Ffat + epsilon2,
age1 = Fage + e11, bmi1 = Fbmi + e12,
fat1 = Ffat + e13,
cholest1 = Fcholest + e14, diastol1 = Fdiastol + e15,
age2 = Fage + e21,bmi2 = Fbmi + e22,
fat2 = Ffat + e23,
cholest2 = Fcholest + e24, diastol2 = Fdiastol + e25;
variance /* Variances of exogenous vars will be
called v-something. __ means fill in the numbers. */
Fage Fbmi Ffat epsilon1 epsilon2 e11-e15 e21-e25 = 15 * v__ ;
cov /* Covariances: If not mentioned, it's zero. */
Fage Ffat Fbmi = 3 * phi__ , epsilon1 epsilon2 = psi12 ,
e11-e15 = 10 * omega1__ , e21-e25 = 10 * omega2__ ;
/* If you don't count the variances and covariances you get a warning.
It's better to count them. */
bounds 0.0 < v01-v15; /* Variances are positive */
/* Wald test of H0: beta12=beta22=0:
Testing the functions f1 = 0 and f2 = 0. Non-linear
functions of the parameters are okay. */
simtests BMI = [f1 f2];
f1 = beta12;
f2 = beta22;