/*

rate       rate per 1000
pop        population
injured    number injured
int_area   indicator of intervention area

*/

*  ------- for SAS ---------   ;

* -- see below for Stata ----- ;


options ls=95 ps=55;

DATA a;
INPUT year rate1000 pop injured int_area;

log_rate = log(injured/pop);

logit = log(injured/(pop-injured));

year_87 = year-1987;

* Make product for interaction term for model c(ii) ;

* Product = ????? *  ????? ; * <-- to be fixed up ;                            ;

*                             and remove the * at
                              beginning of statement ;
LINES;
1983 10.5  3247  34  1
1984 15.0  3200  48  1
1985 10.2  3140  32  1
1986  8.4  3092  26  1
1987 12.4  3056  38  1
1988  7.3  3016  22  1
1989 13.3  3006  40  1
1990  9.1  3072  28  1
1991  9.5  3160  30  1
1983  8.3  4118  34  0
1984  6.9  4046  28  0
1985  8.4  4059  34  0
1986  6.7  4018  27  0
1987  9.5  4018  38  0
1988  4.2  4049  17  0
1989 10.7  4113  44  0
1990 11.0  4179  46  0
1991  7.1  4235  30  0
;
RUN;

PROC PRINT; RUN;

PROC PLOT DATA = a HPERCENT=50;
 PLOT rate1000 * year = int_area
      LOG_RATE * year = int_area / VPOS=25 hpos=30;
RUN;

PROC SORT DATA = a; BY int_area;

PROC REG DATA = a; BY int_area;
 MODEL rate1000 = year; RUN;

PROC REG DATA = a; BY int_area;
 MODEL rate1000 = year_87; RUN;

PROC REG DATA = a;
 MODEL rate1000 =                  ; RUN;

PROC REG DATA = a;
 MODEL rate1000 = year_87          ; RUN;

PROC REG DATA = a;
 MODEL rate1000 = year_87 int_area / CORRB ; RUN;

PROC REG DATA = a;
 MODEL rate1000 = year_87 int_area product  / CORRB ; RUN;


*    ----- Binomial regressions ------ ;


PROC LOGISTIC DATA = a ;
    MODEL injured/pop =                  ; RUN;

PROC LOGISTIC DATA = a ;
    MODEL injured/pop = year_87          ; RUN;

PROC LOGISTIC DATA = a ;
    MODEL injured/pop = year_87 int_area ; RUN;

* TO BE COMPLETED ..  ;

RUN;




* ----------- for Stata  -------------


/*

rate       rate per 1000
pop        population
injured    number injured
int_area   indicator of intervention area

*/

clear

input year rate1000 pop injured int_area
1983 10.5  3247  34  1
1984 15.0  3200  48  1
1985 10.2  3140  32  1
1986  8.4  3092  26  1
1987 12.4  3056  38  1
1988  7.3  3016  22  1
1989 13.3  3006  40  1
1990  9.1  3072  28  1
1991  9.5  3160  30  1
1983  8.3  4118  34  0
1984  6.9  4046  28  0
1985  8.4  4059  34  0
1986  6.7  4018  27  0
1987  9.5  4018  38  0
1988  4.2  4049  17  0
1989 10.7  4113  44  0
1990 11.0  4179  46  0
1991  7.1  4235  30  0
end

gen log_rate = log(injured/pop)

gen logit = log(injured/(pop-injured))

gen year_87 = year-1987

* Make product for interaction term for model c(ii) 

* Product = ????? *  ????? ; * <-- to be fixed up

*                             and remove the * at
*                             beginning of statement 

list

sort int_area

* JH hasn't yet figured out how to overlay two plots
* on same graph.. suggestions welcomed..

by int_area: plot rate1000 year 

by int_area: reg  rate1000 year

by int_area: reg  rate1000 year_87

reg rate1000

reg rate1000 year_87

reg rate1000 year_87 int_area 

reg rate1000 year_87 int_area product
vce, corr


*    ----- Binomial regressions ------ ;


glm injured year_87           , family(binomial pop) link(logit)

glm injured year_87 int_area  , family(binomial pop)  link(logit)

* TO BE COMPLETED ..