options ls=79 ps=50 nocenter;

/* read data (level is score as a factor) */
data dat;
   infile 'damage.dat' firstobs=2;
   input loss score yield;
   level = score;
   logit = log(loss/(1-loss));

/* loss score against yield */
proc glm data=dat;
   class score;
   model yield = loss;
   output out=loss p=pyield r=ryield;
proc plot;
   plot yield*loss=score pyield*loss='*' / overlay;

proc sort; by score;
proc means noprint; by score;
   var yield;
   output out=means mean=myield;
proc glm data=dat;
   class score;
   model yield = score;
   output out=score p=pyield r=ryield;
proc sort; by score;
data merge; merge score means; by score;
proc plot;
   plot yield*score myield*score='-' pyield*score='*' / overlay;

/* test for linearity (\myexpage{plreg}) */
proc glm data=dat;
   class level;
   model yield = score level;

/* score error in variables (\myexpage{plerr}; \myfigpage{plerr}) */
proc glm data=dat;
   class score;
   model yield = score logit;
   output out=both p=pyield r=ryield;
proc plot;
   plot yield*logit=score pyield*logit='*' / overlay;

proc glm data=dat;
   class score;
   model yield = logit;
   output out=logit p=pyield r=ryield;
proc plot;
   plot yield*logit=score pyield*logit='*' / overlay;