/*
Filename:   Pueblos.sas
Purpose:    SAS code for regression analysis of data in Dohm 1990.
Last Update: FDN 3.30.02
*/
title1 'Pueblo Room Counts and Population';
goptions reset=global gunit=in
         ftitle="garamond" ftext="garamond" htitle=.5 htext=.3;

data pueblos;
input Pueblo $ MapDate Population Rooms;
if mapdate = 1882 then period='1884 ';
if mapdate =1948 or mapdate=1952 then period='c1950';
if mapdate > 1970 then period ='1970s';
logpop=log(population);
logrooms=log(Rooms);
cards;
Acoma   1948    879 387
Cochiti 1952    444 225
Isleta  1948    1470 804
Jemez   1948    883 459
Laguna  1948    711 189
Nambe   1948    155 34
Picuris 1948    130 59
SanFelipe   1948    784 276
SanIldefonso    1948    170 207
SanIldefonso    1973    413 189
SanJuan 1948    768 176
Sandia  1948    139 80
SantaAna    1948    288 136
SantaAna    1975    498 152
SantaClara  1948    573 144
SantoDomingo    1948    1106    500
Shipaulovi  1882    113 131
Shongopavi  1882    216 253
Sichomovi   1882    104 105
Taos    1948    907 543
Taos    1973    1463    627
Tesuque 1948    160 116
TewaVillage 1882    175 158
Walpi   1882    270 363
Zia 1948    267 126
;

symbol1 color=red value=dot height=.2 ;
symbol2 color=blue value=dot height=.2 ;
symbol3 color=green value=dot height=.2 ;


proc gplot;
    plot logpop*logrooms =period;
run;

proc glm;
class period;
model logpop= period logrooms /solution;
    output out=regresults
         p=predicted
         r=residual
         uclm=MeanUpperCL
         lclm=MeanLowerCL;


proc sort data=regresults;
    by logrooms;

/*plot the results -- including the residuals! */
proc gplot data=regresults;
    plot residual*predicted = period;
run;

proc univariate data=regresults ;
    var residual;
    qqplot residual /normal (mu=est sigma=est);

goptions reset=global;
symbol1 color=red value=dot height=.2;
symbol2 color=black interpol=join height=2;
symbol3 color=black interpol=join height=2;
symbol4 color=black interpol=join height=2;

proc sort;
by period logrooms;
proc gplot data=regresults uniform;
by period;
    plot    logpop*logrooms
            predicted*logrooms
            MeanUpperCL*logrooms
            MeanLowerCL*logrooms/overlay ;
run;