Homework #3 Solutions STAT 506 Spring 2013

1. We start with an analysis data on rats where the experimental units are 12 thirsty albino rats who are trained to press a lever to get water prior to the experiment. Their pre-experiment pressing rate is recorded as low (1), medium (2), or high (3). They are then injected with one of four levels of a drug where 0 is a control saline solution, the other values are mg per kg of the rat's weight. This was a cross-over design replicated twice, so each rat has 8 measurements of postRate (number of lever presses per second), two at each of four drug levels, and the order of treatments was randomized for each rat. Time intervals between treatments were long enough to remove any carryover effects.

   1. First we look at a plot separating the post rate response by preRate and dose.
      

I see that prerate has a strong effect, dose only kicks in at the highest level, and there is little or no interaction.

2. Write out a model for these data using preRate as a three-level factor and drug as a four-level factor. Include distributions for all random components (assuming normality throughout).

y_ijk = β ₀ + β_i + τ _j + b _k + ε _ijk where k is ratID number from 1 to 12, i is preRate from 1 to 3, and j is dose from 1 to 4.

b_k ˜ iid N(0, σ²_rat), for each rat,
ε_ijk ˜ iid N(0, σ²)

1. Skip: Fit the above model to the data using PROC GLM in SAS using a repeated statement for ratID. Explain the results.

2. Fit the above model to the data using PROC MIXED in SAS using a random term for each rat. Explain results and compare to those just above.

By fitting a random term for each rat ID, we create a compound symmetric variance-covariance structure. Observations within each rat are correlated, but between rats are uncorrelated. SAS estimates the correlation rat-to-rat variance as 0.005146, compared to underlying variance estimated as 0.001393. Correlation is the former over the sum, 0.787. Prerats 1 (low) and 2 (medium) have negative estimates because they are quite a bit lower than the mean for prerate 3 (high) ( -0.45, SE = .037, and -0.21, SE = 0.37 respectively). The effect of dose shows up because we have contrasts between the high dose of 1.8 and control (0.27, se=0.01), low dose (0.28, 0.01) and dose = 1 (0.28, se = 0.01).

3. These results do agree with those of R.

2. We then examined the Sitka data (from the MASS library in R) on the growth of 79 sitka spruce trees.

   1. Here's a plot of the data, and the code I used to create it.

PROC SGpanel data= sitka;

panelBy treat;

series x= time y = size /group = tree;

run;

2. Use PROC GLM to fit a quadratic model across all the data. Update the model adding treatment effects which allow the intercept, slope, or quadratic coefficients to depend on treatment.

I did this in Proc Mixed to get AIC which is 818

3. Using PROC MIXED

   1. Add AR1 correlation structure (within a tree).

AIC drops to -54.9

2. Add compound symmetric correlation (within a tree).

AIC = 94.1

3. Add symmetric correlation (within a tree).

AIC = -74.1Which surprises me, as all covariances are pretty close to the estimate of CS which is 0.37

Compare the four models using AIC.

3. Reduce the model one term at a time until all terms have small p-values in the t-tests.

   I took out treat*timesq because the pvalue was 0.19. Then all terms had small p-values.

4. Plot the residuals versus fitted and normal quantile plots. Discuss any problems.

6. How does the ozone treatment affect growth of these trees? Does SAS give results just like those of R?

As we saw with R, there is a quadratic pattern to growth and the slope for size depends on the treatment.

3. In a study of soil properties, samples were taken on a 10 point by 25 point grid. Download the data here We'll work with two variables: response Ca (calcium concentration) and predictor pH (low numbers are acidic, high numbers basic, 7 is neither).

   1. Make a scatterplot of the two variables and fit a model for Ca based on pH. (Choose the form of the model based on the scatterplot.) Print the estimated coefficients and discuss the relationship.

Still looks linear.

2. Fit the five forms of spatial correlation available in the nlme library. Compare them with each other and with the original model. Do any of the spatial correlation fits improve AIC by more than 2 units? Which is the best of the 5 fits?

All are way better than the original. With SAS, Exponential is best by a whisker, whereas in R it was spherical, but the two were close in each case.

3. Compare results with those we got last semester in HW6 using R.

We have to add the LOCAL option to get a nugget. Then the SAS output agrees well with R.

 /******************  Rats repeated measure  *****/
PROC Import out = rats datafile="~/Projects/stat506/HW3/drugResponse.csv" 
  dbms=dlm    replace;
  delimiter=',' ;
  getnames=yes;
run;
*PROC means;
*run;
PROC sgpanel data = rats;
   panelby dose / columns = 4;
   vbox PostRate /category= PreRate;
run;
PROC sgpanel data = rats;
   panelby prerate / columns = 3;
   vbox PostRate /category= dose;
run;

PROC mixed data = rats;
   class dose prerate;
   model postrate = prerate dose/ solution ss1;
   random intercept/ subject =  ratID;
   run;

/********************* Tree growth *****/

PROC Import out = Sitka datafile="~/Projects/stat506/HW3/sitka.csv"

dbms=dlm replace;

delimiter=',' ;

getnames=yes;

run;

PROC means;

run;

PROC SGpanel data= sitka;

panelBy treat;

series x= time y = size /group = tree;

run;

data sitka;

set sitka;

timeSq = (time - 202)**2;

run;

Proc Mixed data = sitka;

class treat;

model size = treat|time treat|timeSq;

* random intercept / subject = tree type = ;

run;

Proc Mixed data = sitka;

class treat tree;

model size = treat|time treat|timeSq / htype=1;

repeated / subject = tree type = cs;

run;

PROC Mixed data = sitka;

class treat tree;

model size = treat|time treat|timeSq / htype=1;

repeated /subject = tree type = un;

run;

proc Mixed data = sitka;

class treat tree;

model size = treat|time treat|timeSq / htype=1;

repeated /subject = tree type = ar(1) ;

run;

PROC Mixed data = sitka;

class treat tree;

model size = treat|time timeSq / htype=1;

repeated /subject = tree type = un;

run;

17:29 Sunday, January 27, 2013

No correlation:

Compound symmetric

Autoregressive:

Type 1 Tests of Fixed Effects
Effect	Num DF	Den DF	F Value	Pr > F
treat	1	77	1.67	0.1996
Time	1	77	1275.25	<.0001
Time*treat	1	77	15.27	0.0002
timeSq	1	77	273.85	<.0001

/******************   Soils  spatial ********************************/
PROC Import out = soils datafile="~/Projects/stat506/HW3/soils.dat" 
  dbms=dlm    replace;
  delimiter=' ' ;
  getnames=yes;
run;
PROC means;
run;

PROC SGPLOT data = soils;
  scatter x = Ca  y = pH / markerattrs=(symbol=CircleFilled) 
                  transparency=0.7;
  run;

 Proc mixed data = soils;
   model pH = Ca ;
  run;
    
 Proc mixed data = soils;
  model pH = Ca ;
  repeated /local   subject = intercept type = sp(exp) (column row);
  run;
  
Proc mixed data = soils;
  model pH = Ca ;
  repeated / local  subject = intercept type = sp(lin) (column row);
  run;

  
 Proc mixed data = soils;
  model pH = Ca ;
  repeated / local  subject = intercept type = sp(gau) (column row);
  run;
  
  
 Proc mixed data = soils;
  model pH = Ca ;
  repeated /local   subject = intercept type = sp(sph) (column row);
  run;
  

No Covariance:

Exponential

Linear

Spherical