%include data(fuel);

title 'Effects of Errors-in-predictors';

data fuel;
	set fuel;
	label fuel='fuel/person'
		;
	popdens = log10( pop / area );
	
proc means data=fuel mean std;
	var fuel popdens;

%let std = 0.6059;   *-- true standard deviation of posdens;

*-- create some new predictors with varying amounts of error;
data fuel;
	set fuel;
	popden1 = popdens +  .5 * &std * normal(1234524241);
	popden2 = popdens +   1 * &std * normal(1234524241);
	popden3 = popdens + 1.5 * &std * normal(1234524241);
	popden4 = popdens +   2 * &std * normal(1234524241);

proc corr nosimple noprob;
	var fuel popdens popden1-popden4;

proc plot hpercent=50 vpercent=50;
	plot fuel * (popdens popden1 popden2 popden3);

proc reg outest=coeff noprint;
	m0:  model fuel = popdens;
	m1:  model fuel = popden1;
	m2:  model fuel = popden2;
	m3:  model fuel = popden3;
	m4:  model fuel = popden4;

*-- put the coefficients in one column;
data coeff;
	set coeff;
	select (_model_);
		when('m1') popdens = popden1;
		when('m2') popdens = popden2;
		when('m3') popdens = popden3;
		when('m4') popdens = popden4;
		otherwise;
	end;
	drop _type_ fuel popden1-popden4;

proc print data=coeff;
	title2 'Model summary statistics';

*-- Approximate the Intraclass correlation by  var(x) / (var(x)+var(e));
data coeff;
	set coeff;
	select (_model_);
		when('m1') rsq = (&std**2) / (&std**2 + 0.5**2);
		when('m2') rsq = (&std**2) / (&std**2 + 1.0**2);
		when('m3') rsq = (&std**2) / (&std**2 + 1.5**2);
		when('m4') rsq = (&std**2) / (&std**2 + 2.0**2);
		otherwise  rsq = 1;
	end;
	b = popdens / sqrt(rsq);
	run;
proc print;
	title2 'Coefficients adjusted for unreliability';


*include goptions;
goptions htext=2 htitle=2.5;
title 'Original data';
axis1 order=(-4 to 2);
%ellipses(data=fuel,y=fuel, x=popdens, interp=rl, haxis=axis1);
%gdispla(OFF);
title 'Error in X: 0.5 std';
%ellipses(data=fuel,y=fuel, x=popden1, interp=rl, haxis=axis1);
title 'Error in X: 1.0 std';
%ellipses(data=fuel,y=fuel, x=popden2, interp=rl, haxis=axis1);
title 'Error in X: 1.5 std';
%ellipses(data=fuel,y=fuel, x=popden3, interp=rl, haxis=axis1);
title 'Error in X: 2.0 std';
%ellipses(data=fuel,y=fuel, x=popden4, interp=rl, haxis=axis1);
%gdispla(ON);
%panels(rows=2, cols=2, order=down, first=2);