In the Week 9 live workshop we will work through fitting a generalized linear model (in particular, logistic regression) on the Iberian weather data.

1      Setting up

In addition to these instructions you will need the following, all of which can be downloaded from the ‘Week 9’ tab of the Moodle page:

•      The lecture materials.

•      The SAS programs RainLogistic.sas.

•      The data files RainTemp.dat and stations.csv. The RainTemp.dat file is the same that you used in Week 8: if you already have a copy in the same folder as                    RainLogistic.sas then you don’t need to download it again.

2      PROG REG and PROG GLM

In the self-study session you worked through an example of linear regression as well as an ANOVA model. Spend the first 10-15 minutes of the live workshop discussing within your  groups.

A few specific questions that you might want to discuss:

•      Were you able to fill-in the blanks of the sheep energy script successfully?

•      Looking at the coefficients and the p-values in the linear and quadratic models for the   sheep energy example, what differences and similarities do you observe? Can you think why this might be?

•      In step 3 of the battery life exercise (page 9 of the self-study notes), what are the equations of the lines on the interaction plot and why are they parallel?

•      In step 4 of the battery life exercise, why do you get the same model fits whether you treat Temperature as a factor, or as a continuous covariate with linear and quadratic   terms included?

As always, feel free to ask questions if there are questions you can’t answer after discussing within your breakout groups.

3      PROC GENMOD for generalized linear models: analysis of

Iberian weather data

For the final exercise in this week’s workshop, we’ll return to the Iberian weather data from Week

8. The program is RainLogistic.sas: open it and work through it. The purpose is to model the association between the number of ‘wet’ days in each year and the temperature. One possible     application of this could be to answer the question: if the temperature of this region increases in the future, what might happen to the frequency of rainfall events?

Taking the number of wet days as a response variable, a binomial distribution might be suitable  in the first instance (you could consider each of the 365 days as an independent trial, with some  success probability related to the annual mean temperature). Following from the material on       generalized linear models in Week 6, you might therefore be tempted to try a logistic regression model with temperature as a covariate.

The main SAS procedure for generalized linear models is PROC GENMOD. Having seen both PROC  REG and PROG GLM though, most of it should feel fairly familiar to you: the use of ODS graphics   to obtain nice graphical output, the use of CLASS statements to define factor covariates (don’t     worry, there’s nothing as difficult as the battery example here!) and the way that model formulae work. PROC GENMOD also produces measures like AIC (see Week 6) which can be used to               compare models very easily. You should note the caveats on this from Week 6, however (see        Section 2.1 of the notes from that workshop).

The main new feature for PROC GENMOD is the need to define a distribution and a link function in the MODEL statement, for example:

This fits a model in which:

•      The response variable is the number of wet days (NWet) as a proportion of the total number of days (NDaysR — this is the number of days in a year for which the rainfall value is non-missing);

•      The covariate is temperature (Temp);

•      The response distribution is binomial (DIST=bin) with a logit link (LINK=logit).

To see what other distributions are available and how to define them, look at the help for PROC GENMOD. Apart from the use of PROC GENMOD, there are two parts of the RainLogistic.sas

program that deserve special attention here. They are indicated by “!!! SEE WORKSHOP NOTES !!!” as usual. They are as follows:

•      After reading the daily data, the program creates an annual data set using a MEANS step:

You have seen this kind of thing before (Week 8 notes, Section 6.4). This particular example is slightly more advanced though. Note the following:

–      According to the CLASS statement, the output data set will contain a row for each Site:Year combination — i.e. annual values at each site separately.

–      Each row of the output data set will contain variables Nwet (the number of wet

days), Temp (the mean temperature), NDaysR (the number of non-missing daily    rainfall observations) and NDaysT (the number of non-missing daily temperature observations). These variables are all defined on the second line of the OUTPUT  statement above. Notice how the OUTPUT statement allows us to produce            different statistics for different variables — the totalnumber of wet days but the   meantemperature, for example.

–      The OUTPUT statement also allows us to select a subset of the rows to include in

the output data set. This is done via the WHERE secondary option. You saw the use of the _TYPE_ argument in Week 8, so you should know roughly what it does. If    you’re not sure, remove it from the code and see what happens (then make sure   you put it back!). The NDaysR>360 & NDaysT>360 condition specifies that the

output data set should be restricted to cases where there were at least 360 non-   missing daily values for both rainfall and temperature. The reason for this is that if there were substantial numbers of missing daily values at some site in a particular year, they may all be in one season and hence bias the results.

•      After fitting Model 4 (a quasibinomial generalized linear model), the program imports

some additional information on the stations — their geographical co-ordinates (latitude   in degrees north of the equator and longitude in degrees east — the longitude values are all negative, therefore) and altitude in hundreds of metres. This additional information is  then combined with the existing data using PROC MERGE — it’s often useful to be able to  do this kind of thing, although with SAS you need to ensure that the data sets are all        sorted in the right order before merging them. Apart from the use of PROC MERGE, the     main point to note is that the station information is stored in a CSV file, instead of the      usual space-delimited file (this is because the file also contains station names, which         themselves contain spaces). To specify that the values are separated by commas rather     than spaces, the INFILE statement specifies the DLMSTR="," (‘delimiter string is a            comma’) option.

If you look at the CSV file,1you’ll see that in fact it contains much more information about each station than just the geographical coordinates and altitude. If you want to explore this in more  detail — e.g. by importing more station attributes to see if these can be related to the numbers of wet days at each station — you’re welcome to do so!

This exercise should take you about 45 minutes. When you have finished, as well as having a bit more experience with GLMs you should be able to answer these questions:

•      If temperatures in the region increase in the future, would you expect any change in the frequency of rainfall events?