Assignment 8

Due Date

Friday, March 16, 2012

Data Source

The data set for this assignment is burn.csv, a comma-delimited text file.

Overview

The data come from a study of the role that prescribed burning plays as a land management tool in Onslow Bight, a region of the North Carolina (NC) coastal plain where a multiagency partnership has been established for the conservation of the longleaf pine ecosystem. The objective was to determine the site and landscape characteristics that best explain the placement of prescribed burns as a way to assess whether prescribed burning is being carried out in areas that need restoration.

The authors compiled GIS data delineating the locations of prescribed burns conducted between 1989 and 2007 from the six Onslow Bight land-management agencies that manage their land via prescribed burning: Camp LeJeune, Croatan National Forest, Cherry Point, Cedar Island NWR, NC WRC, and TNC. Soil survey data were used to select only sites that are predominantly on soils with relatively low organic matter content—the soils that support longleaf pines. This identified sites that currently support longleaf pine communities or that could be restored to longleaf and yielded 697 sites covering a total of 52,226 hectares.

Records of prescribed burning were in the form of GIS polygon data delineating burned areas from 1989 to 2007. Sites were included in the study only if they provided at least 15 years of burn history data. Analysis was then restricted to the years 2004 to 2007. For each year during the period 2004 to 2007 the authors labeled a site as “burned” for that year if 50% or more was burned. Thus each of the 697 sites in the study contributed four binary observations.

The full history of prescribed burning, beginning in 1989, was used to create the set of predictor variables that were thought to influence whether a site was burned. The authors considered a number of different models. We will focus on a model that included only what the author's termed non-ecological variables. These include:

1. burn50 is a binary variable indicating whether more than 50% of the site was burned in a given year. This is the response variable of the analysis.
2. year records the year of the measurement.
3. PolyFID is a categorical variable that identifies the site at which the four annual burn measurements were obtained.
4. Hectares is the area of the site.
5. zmin_drd is the distance from the nearest road.
6. zmin_ddv is the distance from the nearest development.
7. zmin_dr is the distance from the most recent burn.
8. prop_timb is the proportion of the site in managed timber.
9. Prop_no_b is the proportion of the site that has never been burned.
10. rx_cat3 is a categorical variable with four categories that denote the time since the last burn: 1= 1 year, 2 = 2-3 years, 3 = 4-5 years, and 4 = greater than 5 years.

Some of these variables varied by year at each site (e.g., time since last burn and distance from the nearest burn) while others were treated as constant at a site during the four years covered by the study.

The data in this study consist of repeated measurements on the same site, as well as separate measurements on different sites. Thus there is the potential for temporal correlation among observations coming from the same site. A site’s burn status in a given year can affect the decision to burn in following years. In addition, sites are aggregated in space, leading to potential spatial correlation if the decision to burn one site affects the decision to burn a nearby site. In this homework assignment we will address only the temporal correlation in these data.

Questions

1. Fit a model with the binary variable burn50 as the response. As predictors use Hectares, rx_cat3 (as a factor), zmin_drd, prop_timb, Prop_no_b, zmin_dr, zmin_ddv, and the interaction of zmin_ddv and rx_cat3. For this preliminary model ignore the issue of temporal correlation.
2. Use GEE to include a temporal correlation model for the response. Use the same set of predictors as in Question 1 but add the following correlation structures:
   1. independence
   2. exchangeable
   3. unstructured
   4. AR(1)
   5. stationary M-dependent of order 1
   6. non-stationary M-dependent of order 1
3. Use appropriate statistical criteria to determine which of the six correlation models is best for these data.
4. A common rule-of-thumb for management of longleaf pine ecosystems is to burn every three years. Assuming that at least some of the managers were following this rule of thumb explain why it gives additional support for the correlation model you selected in Question 3.
5. Fit a generalized linear mixed effects model with random intercepts to these data using the same set of regressors that you used in Question 1 (and Question 2).
6. Assemble in a table the estimated coefficients for the generalized linear model, the best generalized estimating equation model, and the fixed effect estimates from the generalized linear mixed effects model. Two of the sets of estimates will be very similar. Explain why this probably happened and in particular why the third set of estimates is different from these two.
7. Examine the estimates of the elements of the working correlation matrix for the model you chose in Question 3. Determine which if any of the entries fall outside the Prentice bounds for correlations between binary observations.
8. Replace the correlations in the correlation matrix that fall outside the Prentice bounds with values that lie just inside the Prentice bounds. Use two-digit precision for this, i.e., the new correlation should be recorded to only two significant digits yet still fall within the bounds. Leave the remaining correlations alone.
9. Refit the GEE model using your new correlation matrix treating the correlation matrix as fixed. Report the summary table of the model.

Hints

• General: Some of the GEE questions will need to be answered using the gee package while others will need to be answered using the geepack package.
• Question 7: A random intercepts model induces a correlation in the response variable. What is the nature of that correlation and why does knowing that help you answer this question?
• Questions 8: This material was not covered in the formal lab session but the information for doing this is described in the notes for lecture 23. The situation you have here is much simpler than it was in lecture 23 because you have no missing data and each site provides data for all four years. Thus you can use the raw response and the fitted values from the model directly and can skip over the steps where missing values were introduced to fill out the groups. You have six different correlation values to examine. You will need to construct a two-column matrix that contains the six different ways of pairing up the four years. You can either use the for loop code from the notes or construct the matrix by hand.
• Question 9: The details for doing this are explained at the end of the lecture 23 notes. It was not carried out there because the geeglm function could not fit the model without crashing R. You won't have that problem here. You will need to specify corstr="fixed" and create a zcor object that contains the six unique non-diagonal entries of the correlation matrix (the lower triangle treated in column order) as a vector. The entire vector then must be repeated as many times as there are unique groups in the data set. As a result the zcor object should have length that is six times the number of rows in the data set. See the help screen of the R geese function for an example.

Course Home Page