Lecture 27 (lab 7)—February 28, 2006

What was covered?

• R topics
  • graphics parameter options mar and las
  • sorting data frames by a specific column using order
  • residual plots with cr.plots
  • controlling glm and glm.nb iterations with the control argument
• checking model form with residuals
• offsets and covariates in generalized linear models
• checking model residuals for spatial correlation

R functions and commands demonstrated

• as.matrix converts a data frame into a matrix
• cr.plots (from car library) generates component plus residual (partial residual) plots
• glm.control a function for setting options on the numerical estimation routine of glm
• order returns a vector of the current indices of elements of a vector but arranged as if the elements themselves were sorted
• qres.nbinom (from statsci.org) creates randomized quantile residuals for a negative binomial regression model
• unique returns the unique elements in a vector (removing duplicates)

R function options

• control= (argument to glm or glm.nb) for changing default options of the numerical estimation method. Used in conjunction with glm.control
• las= (argument to par) can be used to specify the orientation of the labels for tick marks. The default setting is las=0 which orients the text to be parallel to the axis. Specifying las=1 forces all labels to be horizontal.
• mar= (argument to par) specifies the settings for the four margins of a graph. The four settings are entered in the following order: bottom, left, top, and right. The default settings are 5.1 4.1 4.1 2.1 measured in "line" units
• maxit= (argument to glm.control) specifies the maximum number of iterations to be used in the numerical estimation algorithm
• type=deviance (argument to residuals) is used to extract deviance residuals from a generalized linear model object
• type=pearson (argument to residuals) is used to extract Pearson residuals from a generalized linear model object

R packages used

• car for the cr.plots and qq.plot functions.
• MASS for the function glm.nb The MASS package implements procedures described in Venables and Ripley (2002)
• date for the function as.date
• qres for the qres.nbinom function. Not a package but an online repository of code for computing randomized quantile residuals. http://www.statsci.org/s/qres.html
• residual.NB.MoranI and plot.NB.MoranI are two functions I've written for carrying out a randomization test of Moran's I for spatial autocorrelation of the residuals from a negative binomial regression. http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/Moran residual.doc
• get.acf is a function I've written to calculate the autocorrelation function for a variable based on a list of actual lags. See http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/ACF function.txt

Overview

• We revisit the coral reef data set that we considered during the first few weeks of this course. At that time we investigated the relationship between disease prevalence in various reefs (as measured by the number of observed diseased colonies, a count) and a temperature stress metric that the researchers called WSSTA. We drew the following conclusions from these data.
  • The mean and variance of disease prevalence seem to be related. If we stratify the observations by the level of WSSTA we find that the variance is well described by a quadratic model in the mean. As we've learned since this is consistent with assuming a negative binomial distribution for the random error component.
  • Preliminary loess plots suggested that disease prevalence monotonically increases with WSSTA up until a point, after which it decreases. This suggests a quadratic model rather than a simple linear model may be more appropriate.
  • The original data exhibit a spatial and temporal structure. Observations are spatially clustered and the observations in these clusters were sampled at roughly the same time.
• In this session we attempt to deal with these issues and others by fitting appropriate models and then assessing the results via residual analysis

Fitting the model

> #obtain data
> corals<- read.table('http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/year5coral.txt', header=TRUE, sep=',')

• Although preliminary evidence suggests a negative binomial regression model might be a sensible starting point, I elect to begin by fitting a Poisson regression model to establish a benchmark for comparison.

> model0<-glm(PREV_1~WSSTA,data=corals,family=poisson)
> summary(model0)

Call:
glm(formula = PREV_1 ~ WSSTA, family = poisson, data = corals)

Deviance Residuals:
    Min     1Q Median    3Q    Max
-17.790 -5.032 -2.329 2.510 16.088

Coefficients:
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 2.198401   0.051729   42.50   <2e-16 ***
WSSTA       0.268846   0.006051   44.43   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 4660.4 on 46 degrees of freedom
Residual deviance: 2716.8 on 45 degrees of freedom
AIC: 2914.9

Number of Fisher Scoring iterations: 6

• The Poisson regression results indicated that WSSTA has a significant relationship to prevalence. The estimated coefficient of WSSTA is 0.268846. Since our link function is a log, this coefficient reflects changes in the log mean prevalence. The reported model is: log μ = 2.198 + 0.269*WSSTA. To understand the effect of WSSTA on the original scale we need to exponentiate this equation.

Thus we see that each one unit increase in WSSTA multiplies the predicted mean prevalence by 1.3

• Unfortunately we can also see from the output that this model does not fit. For a Poisson model the residual deviance can be used as a goodness of fit statistic. The reported residual deviance is 2716.8 on 45 degrees of freedom. For a good-fitting model the ratio of the residual deviance to its degrees of freedom should be approximately 1. Since the ratio here is far greater than 1, we conclude the data are overdispersed. This is additional evidence that a negative binomial error distribution would be a more appropriate choice here.

> library(MASS)
> model1<-glm.nb(PREV_1~WSSTA,data=corals)
> summary(model1)

Call:
glm.nb(formula = PREV_1 ~ WSSTA, data = corals, init.theta = 0.521663980407755,
link = log)

Deviance Residuals:
    Min      1Q  Median     3Q    Max
-2.2950 -1.1682 -0.3903 0.2552 1.5750

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.13280    0.37132   5.744 9.26e-09 ***
WSSTA        0.28231    0.06564   4.301 1.70e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.5217) family taken to be 1)

    Null deviance: 75.547 on 46 degrees of freedom
Residual deviance: 55.937 on 45 degrees of freedom
AIC: 413.98

Number of Fisher Scoring iterations: 1

Correlation of Coefficients:
(Intercept)
WSSTA -0.84

Theta: 0.522
Std. Err.: 0.102

2 x log-likelihood: -407.981

• Observe that the AIC has decreased tremendously, from 2914.9 in the Poisson model to 413.98 in the negative binomial model. Interestingly there has been little change in the parameter estimates. Thus our interpretation of the effect of WSSTA on the mean prevalence is roughly the same as it was for Poisson model.
• Since preliminary graphical analysis suggested a quadratic model might be appropriate, I try fitting that next. Instead of creating the squared term explicitly in the data frame I create it on the fly by sandwiching an expression for it inside the I (identity) function and including it this way in the model formula. Recall that use of the I function allows arithmetic to be carried out inside a model formula. The effect here is to square the variable WSSTA before the model is estimated.

> model2<-glm.nb(PREV_1~WSSTA+I(WSSTA^2),data=corals)
> summary(model2)

Call:
glm.nb(formula = PREV_1 ~ WSSTA + I(WSSTA^2), data = corals,
init.theta = 0.526918648778831, link = log)

Deviance Residuals:
    Min      1Q  Median     3Q    Max
-2.3293 -1.1626 -0.4022 0.4045 1.7069

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.74721    0.53704   3.253 0.00114 **
WSSTA        0.45633    0.21039   2.169 0.03008 *
I(WSSTA^2)  -0.01400    0.01742  -0.804 0.42151
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.5269) family taken to be 1)

    Null deviance: 76.239 on 46 degrees of freedom
Residual deviance: 55.842 on 44 degrees of freedom
AIC: 415.39

Number of Fisher Scoring iterations: 1

Correlation of Coefficients:
          (Intercept) WSSTA
WSSTA     -0.87
I(WSSTA^2) 0.72       -0.95

Theta: 0.527
Std. Err.: 0.103

2 x log-likelihood: -407.393

• Surprisingly and contrary to expectations the quadratic term is not statistically significant and the AIC has actually increased. Thus the statistical evidence we've obtained belies our earlier graphical evidence on the need for a quadratic term.

Checking to see if observations are equivalent

Fig. 1 Coral Cover Distribution

• When count data record the occurrence of some event over space, the question arises as to whether the observations are comparable or whether some ad hoc adjustment needs to be made. For example, if the reefs being examined were of different areas we would expect this to affect the number of events we would see and we would want to compensate for this. For the reefs in this data set standard-sized transects were used so that the amount of area examined did not differ between the reefs.
• Different reefs can be patchy in the number of colonies present. Included in the corals data set is a variable CORAL_COVE that records the percentage of the reef that is covered by coral. Let's examine this variable to see if it shows appreciable variability across the various reefs.
• The following code produces the plot shown in Fig. 1. Observe that coral cover does vary widely among the various reefs although less so among reefs drawn from the same geographic sector. Clearly we need to incorporate these differences in coral cover in our model.

>#sort the reefs by geographic sector
> sortedcorals<- corals[order(corals$SECTOR),]
>#obtain a list of SECTORs
> sector.list<- unique(sortedcorals$SECTOR)
>#save the current margin settings
> par("mar")->oldmar
>#allow for a wider left margin (2nd entry)
> par(mar=c(5.1,6.1,2.1,2.1))
>#change axis label orientation to horizontal
> par(las=1)
>#plot sorted cover values in order
> plot(sortedcorals$CORAL_COVE,1:47, axes=FALSE, xlab='% coral cover', ylab='')
>#add labels that identify the reefs
> axis(2, at=1:47, labels=as.character(sortedcorals$REEF_NAME), cex.axis=.5)
> axis(1)
> box()
>#add horizontal lines to separate observations
> abline(h=1:47, col='grey80', lty=4)
>#use loop to color the points by the geographic sector they occupy
> for(i in 1:6) {
   cur.dat<-sortedcorals[sortedcorals$SECTOR==sector.list[i],]
   cur.num<-(1:47)[sortedcorals$SECTOR==sector.list[i]]
   points(cur.dat$CORAL_COVE, cur.num, pch=16, col=i)
}
> legend(65,47, sector.list, col=1:6, cex=rep(.7,6), bty='n', pch=rep(16,6))
> par(mar=oldmar)
> par(las=0)

Explanation of the code used to produce the graph

• The first line sorts the reefs by the variable SECTOR. The preferred way to do this is with the order function. This function doesn't actually sort anything but instead returns a list of coordinates for the observations as shown below.

  > order(corals$SECTOR)
   [1]  7 17 19 20 30 32 34 39 43 44  9 26 38 46 10 14 28 29 31 33 36 37
  [23] 47  5  6 12 16 18 23 45 11 13 15 21 25 35 40 42  1  2  3  4  8 22
  [45] 24 27 41

• The way to interpret this is as follows. If the observations were rearranged so that 7th observation was followed by the 17th which was followed by the 19th, etc. then the observations would be ordered by SECTOR. We can use this by entering these numbers as the row indices for the data frame corals. What will be returned are the rows of the data frame but now arranged in this new order. The effect is to sort the corals data set by SECTOR.

> sortedcorals<-corals[order(corals$SECTOR),]

• The second line pulls out the unique values of the variable SECTOR into a variable called sector.list. The unique function removes all the duplicate values from a vector. So when we examine this variable we see it only has six different values.

> sector.list<- unique(sortedcorals$SECTOR)
> sector.list
[1] CA CB CL SW TO WH
Levels: CA CB CL SW TO WH

• The next two lines change some of the default graphics settings. Because I plan to include the reef names as the y-axis tick mark labels, I need to expand the left margin of the graph. This is done by specifying new values to mar (for margins) using the par function. The default settings of mar are shown below.

> par("mar")
[1] 5.1 4.1 4.1 2.1

• What's displayed are the current margin settings starting at the bottom of the graph and moving clockwise. The values denote the number of lines from the edge of the graph. The second value corresponds to the left side of the graph and is the one that needs to be changed. To do so just give mar a new vector of values with a larger second value. I also choose to reduce the margin at the top. In order that the default values for the margins can be restored later I first copy the current settings to a vector a I call oldmar.

> par("mar")->oldmar
> par(mar=c(5.1,6.1,2.1,2.1))

• By default axis tick mark labels are printed parallel to the axis they label. This corresponds to the default value las=0. I change this to las=1 which causes all axis tick mark labels to be printed horizontally.

> par(las=1)

• There are lots of additional graphics options. The full set is displayed in the help window for the par function.

• The next line plots the data. By specifying the y-values as 1:47 R just plots the cover values from bottom to top in the order they occur in the variable sortedcorals$CORAL_COVE. I set axes=FALSE because I want to supply my own labels for the tick marks.

> plot(sortedcorals$CORAL_COVE,1:47, axes=FALSE, xlab='% coral cover', ylab='')

• The first axis statement creates tick marks for each of the 47 reefs and prints the reef names at the tick marks along the y-axis. The as.character function is used to ensure that the text of the reef name is displayed, rather than its number. This is necessary because REEF_NAME was by default converted to a factor variable when the data set was read into R with the read.table function.

> axis(2, at=1:47, labels=as.character(sortedcorals$REEF_NAME), cex.axis=.5)
> axis(1)
> box()

• The abline function adds horizontal lines at each reef to visually facilitate associating a plotted point with a specific reef.

> abline(h=1:47, col='grey80', lty=4)

• I use a for loop to color the plotted points so that different colors are used depending on which geographic sector the reef occurs in. Since I previously sorted the reefs by sector, the different colors will cluster together on the plot.
  • In order that multiple lines of code are treated as being part of the same loop, the loop is started with a left-facing curly brace, {, and ends with a right-facing curly brace, }.
  • The first line of the body of the for loop subsets the data set. The variable sector.list contains the names of the different sectors. Specifying sector.list[i] pulls out a sector name depending on the value of i, the index variable of the for loop. Thus the first time through the loop i = 1 and sector.list[i] will be CA. The Boolean expression sortedcorals$SECTOR==sector.list[i] is tested for each observation and evaluates to TRUE or FALSE. It evaluates to TRUE if an observation comes from the sector currently being considered. As a result these are the only observations that get assigned to the data frame cur.dat. So, e.g., when i = 1, cur.dat only contains observations from sector CA.
  • The second line of the loop subsets the numbers 1 to 47 to pull out those reef numbers that correspond to reefs in the currently selected sector.
  • The third line of the loop plots the coral cover values for the current sector. The color selection is from R's default numeric codes for colors. Since i varies from 1 to 6, the colors that are used are the color codes 1 to 6. It would be easy to customize this by creating a list of six colors in a vector, call it color.codes for example, and then replacing col=i with col=color.codes[i].

  > for(i in 1:6) {
     cur.dat<-sortedcorals[sortedcorals$SECTOR==sector.list[i],]
     cur.num<-(1:47)[sortedcorals$SECTOR==sector.list[i]]
     points(cur.dat$CORAL_COVE, cur.num, pch=16, col=i)
     }

• The first line outside of the for loop creates the legend. Each legend component is a vector with six elements because there are six different sectors being displayed.

> legend(65,47, sector.list, col=1:6, cex=rep(.7,6), bty='n', pch=rep(16,6))

• The final two lines reset the graphics parameters back to their default values. The old margin settings were saved in the variable oldmar and are now used to return the margin widths to their defaults.

> par(mar=oldmar)
> par(las=0)

Dealing with observations that are not equivalent

• The graphical display of Fig. 1 suggests that different reefs vary considerably in their coral cover. Intuitively it seems rather obvious that coral cover has the potential of modifying the relationship between disease prevalence and WSSTA. Obviously reefs with more coral cover have more colonies and if there are more colonies there are more colonies that can get sick. Does this pattern actually manifest itself in the data set? Fig. 2 addresses this question.

> plot(corals$CORAL_COVE, corals$PREV_1, xlab='% coral cover', ylab='Prevalence')

• It's clear from the plot that coral reefs with a higher percentage of coral cover tend to have a higher disease prevalence. We need to control for coral cover in our model.
• Lecture 24 addressed the use of control variables in count data regression models. There we discussed two options:
  1. Include the variable as a regressor in the model, log-transformed or not, in which case we call it a covariate.
  2. Include the log of the variable in the model as an offset.
• When a variable is included as an offset, we don't estimate its coefficient but instead hold it fixed at 1. As explained in Lecture 24, the effect is to convert the response variable of the regression model from a count into a rate. This approach can make sense when the rate is readily interpretable as the number of events per unit time or per unit area. With a variable such as we have here, percent coral cover, it's not clear that the notion of a rate makes any sense. As a result, I compare three different ways to include coral cover in the model: log-transformed as a covariate, log-transformed as an offset, and untransformed as a covariate.

> glm.nb(PREV_1~WSSTA+log(CORAL_COVE),data=corals)->model2a
> glm.nb(PREV_1~WSSTA+offset(log(CORAL_COVE)),data=corals)->model2b
> glm.nb(PREV_1~WSSTA+CORAL_COVE,data=corals)->model2c

• If we examine the estimated coefficient of log(CORAL_COVE) in model2a, we see it is far from the fixed value of 1 that is assumed in the offset model, although if we were to construct a Wald confidence interval for it we would find that the interval would just include 1 (details not shown).

> coef(model2a)
(Intercept)     WSSTA log(CORAL_COVE)
  0.2173891 0.2199193       0.6603233

• Using AIC to compare the models we see that the AIC-best model treats coral cover as a covariate that is not log-transformed. The reported AIC is a substantial decrease from the AIC of a model that does not include coral cover (model 1 above, AIC = 415.39).

> AIC(model2a,model2b,model2c)
        df      AIC
model2a  4 407.4004
model2b  3 408.3810
model2c  4 403.3482

Revisiting the structural form of the model

• Model residuals can be used to determine if a predictor has been entered into a model correctly. Since the residuals represent that part of the response that is not explained by the model, a plot of the residuals against predictors already in the model should reveal no additional pattern if the relationship has been correctly specified. From lecture 24 we learned that there are two kinds of residuals that are commonly used as diagnostic tools, Pearson residuals and deviance residuals, of which deviance residuals are usually preferred. Because it is often difficult to detect patterns in residual plots, it is useful to add a lowess curve to the display. Fig. 3 shows both Pearson and deviance residuals with a lowess curve superimposed. As the plot reveals there is a slight hint of a quadratic pattern to the residuals as a function of WSSTA.

> plot(corals$WSSTA, residuals(model2c, type='pearson'), xlab='WSSTA', ylab='Pearson residual')
> lines(lowess(residuals(model2c, type='pearson')~corals$WSSTA), col=2, lty=2)
> plot(corals$WSSTA, residuals(model2c, type='deviance'), xlab='WSSTA', ylab='Deviance residual')
> lines(lowess(residuals(model2c, type='deviance')~corals$WSSTA), col=2, lty=2)

Fig. 3 Model residuals (Pearson and deviance) plotted against WSSTA

• A graphical device that is sometimes more informative than a straight residual plot is something called a component plus residual plot (also called a partial residual plot). A partial residual is a residual in which the linear effect of a specific predictor is added back. The partial residual e_ij is defined as follows.

Here e_i is the ordinary residual for observation i, x_ij is the value of predictor x_j for observation i, b_j is the estimated regression coefficient of predictor x_j, and e_ij is the value of the partial residual for predictor x_j for observation i. These are then plotted against x_j just as in Fig. 3. In our example, x_j is WSSTA.

• The car package contains the function cr.plots that produces component plus residual plots through an interactive menu. It adds a least squares line to the plot that is supposed to be interpreted as the regression surface viewed on edge in the direction of x_j. It also plots a nonparametric smooth that can be contrasted with the least squares line for assessing model specification. For the first plot I accept the default span = 0.5 (Fig. 4). The span corresponds to the fraction of the data set that is used in generating the lowess regression line at a given point.

> library(car)
> cr.plots(model2c)

1: Change span = 0.5
2: WSSTA
3: CORAL_COVE

Selection: 2

1: Change span = 0.5
2: WSSTA
3: CORAL_COVE

• For the second plot I reduce the span to 0.4.

Selection: 1
span: .4

1: Change span = 0.4
2: WSSTA
3: CORAL_COVE

Selection: 2

1: Change span = 0.4
2: WSSTA

          

Fig. 4  Component plus residual plots with different spans: span = 0.5 (left) and span = 0.4 (right)

• This last plot does suggest that perhaps there is some residual nolinearity. I try adding a quadratic term in WSSTA to a model with WSSTA and coral cover.

> model3<-glm.nb(PREV_1 ~ WSSTA + CORAL_COVE + I(WSSTA^2), data = corals)
Warning messages:
1: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
2: alternation limit reached in: glm.nb(PREV_1 ~ WSSTA + CORAL_COVE + I(WSSTA^2), data = corals)

• From the warning message we see that the model failed to converge. We can increase the maximum allowed iterations from its default value of 25. This can be done with the control argument in which we use the glm.control function to change the value of maxit from 25 to 100.

> model3<-glm.nb(PREV_1 ~ WSSTA + CORAL_COVE + I(WSSTA^2), data = corals, control=glm.control(maxit=100))
> summary(model3)

Call:
glm.nb(formula = PREV_1 ~ WSSTA + CORAL_COVE + I(WSSTA^2), data = corals,
control = glm.control(maxit = 100), init.theta = 0.767698858801052,
link = log)

Deviance Residuals:
    Min      1Q  Median     3Q    Max
-2.2579 -0.9526 -0.3898 0.1286 2.5436

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.197456   0.569032  -0.347 0.728589
WSSTA        0.709989   0.183987   3.859 0.000114 ***
CORAL_COVE   0.050201   0.009937   5.052 4.38e-07 ***
I(WSSTA^2)  -0.047068   0.015578  -3.021 0.002516 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.7677) family taken to be 1)

    Null deviance: 107.104 on 46 degrees of freedom
Residual deviance:  54.713 on 43 degrees of freedom
AIC: 398.48

Number of Fisher Scoring iterations: 1

Correlation of Coefficients:
           (Intercept) WSSTA CORAL_COVE
WSSTA      -0.82
CORAL_COVE -0.58        0.23
I(WSSTA^2)  0.75       -0.95 -0.33

Theta: 0.768
Std. Err.: 0.162

2 x log-likelihood: -388.482

• Observe that the quadratic term is statistically significant and that the AIC has decreased. Both measures suggest we should retain the quadratic term. Thus it wasn't until we controlled for coral cover that the need for a quadratic term could be formally demonstrated.

Additional residual analysis

• Because of the repeated integer values of the response (multiple observations with the same response) deviance and Pearson residuals are generally not normally distributed. Worse they often show banded patterns in plots that detract from the display and make it difficult to detect patterns of interest. Dunn and Smyth (1996) proposed a mild transformation of ordinary residuals to yield what they call randomized quantile residuals. Randomized quantile residuals retain the useful diagnostic properties of ordinary residuals, but lack their detracting features.
• Code for calculating randomized quantile residuals for a number of different generalized linear models is available at the web site http://www.statsci.org/s/qres.s. To use their function for negative binomial regression in R we need to modify one line of their code to correctly identify the location of the estimated negative binomial size parameter. Their code for producing quantile residuals from a negative binomial regression is shown below. The line colored red is the line that I changed from what is found at the web address.

qres.nbinom <- function(glm.obj)
{
# Quantile residuals for Negative Binomial glm
# GKS 22 Jun 97
#
y <- glm.obj$y
size <- glm.obj$theta
mu <- fitted(glm.obj)
p <- size/(mu + size)
a <- ifelse(y > 0, pbeta(p, size, pmax(y, 1)), 0)
b <- pbeta(p, size, y + 1)
u <- runif(n = length(y), min = a, max = b)
qnorm(u)
}

• For these data, it turns out there isn't much of a need to compute the randomized quantile residuals. The ordinary residuals don't show any unusual patterns. As an illustration I compare a normal probability plot of the deviance residuals with a normal probability plot of the randomized quantile residuals. I use the qq.plot function from the car package. For these data the plots are almost identical. I have seen other data where the difference is very pronounced.

> qq.plot(residuals(model3,type='deviance'))
> qq.plot(qres.nbinom(model3))

Fig. 5  Normal probability plots of deviance residuals and randomized quantile residuals

Spatial autocorrelation

• There's a spatial component to these data. Reefs don't occur in isolation but occur near other reefs. Furthermore the sample wasn't selected randomly. Sampled reefs occur in clusters both temporally and spatially. Furthermore the values of the variables of interest, disease prevalence and WSSTA also exhibit this same clustering. This was shown in a series of plots done in Lecture 3 (Figs. 13–14, 16–17). Spatial correlation can affect both the interpretation of parameter estimates and the magnitude of their standard errors.
• It's possible though that much of the spatial correlation exhibited by disease prevalence can be explained through its relationship with WSSTA and coral cover. Thus one place to start in addressing spatial correlation is to fit an ordinary regression model assuming independent observations and then examine the residuals of the regression model to see if any lingering spatial correlation can be detected. This is the approach taken here.
• Spatial correlation is often assessed using a correlogram in which the correlation is calculated for observations at progressively greater distances apart. A common measure of spatial correlation is Moran's I, a weighted correlation coefficient used to detect deviations from spatial randomness.
• Because of the extreme patchiness of the data here, treating distance continuously turns out not to be particularly productive. Instead I treat the observations as lattice data (observations occupying fixed locations rather than being sampled from a continuum) and assign observations to neighborhoods defined by their geographic sector. Pairwise weights were constructed such that observations within the same sector shared a common positive weight, while pairs of observations from different sectors received zero weight. These weights were then assembled in a 47 × 47 matrix. I read the weights into R next.

> wtmat<- read.table('http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/wtmat5.txt')
> wtmat5<-as.matrix(wtmat)
> wtmat5[1:7,1:7]
     V1    V2    V3    V4        V5        V6 V7
1 0.000 0.125 0.125 0.125 0.0000000 0.0000000  0
2 0.125 0.000 0.125 0.125 0.0000000 0.0000000  0
3 0.125 0.125 0.000 0.125 0.0000000 0.0000000  0
4 0.125 0.125 0.125 0.000 0.0000000 0.0000000  0
5 0.000 0.000 0.000 0.000 0.0000000 0.1666667  0
6 0.000 0.000 0.000 0.000 0.1666667 0.0000000  0
7 0.000 0.000 0.000 0.000 0.0000000 0.0000000  0

• The as.matrix function turns the data frame wtmat into a matrix wtmat5. The weights used are binary weights so that all observations in the same sector share the same weight. The software used to create the weights divided these binary entries by the number of neighbors an observation has in the neighborhood to yield the fractions shown in the output above. Since 0.125 is one-eighth, we conclude that observations V1–V4 each have eight neighbors (not including itself), V5 and V6 have 6 neighbors, and apparently the neighbors of V7 do not appear among the first six observations listed.
• Let denote the vector of residuals—raw, Pearson, or deviance—from the regression model, then Moran's I is calculated as follows.

Because of the way the weight matrix W is scaled for these data, the first term in the product is just one, so to compute Moran's I we can just calculate the ratio of quadratic forms that appear in the second term.

> #calculate 1'W1
> rep(1,47)%*%wtmat5%*%rep(1,47)
    [,1]
[1,] 47
> #calculate n
> dim(wtmat5)
[1] 47 47

• While formulas for the moments of Moran's I have been derived and under certain conditions Moran's I has an asymptotic normal distribution, the statistical significance of Moran's I is typically assessed with a randomization test. The randomization procedure can proceed under what Cliff and Ord (1981) call assumption R or assumption N.
  • Under assumption R we would randomly permute the residuals, or equivalently, the original observations. The problem with this approach is that if there is any correlation among the regressors this will be disrupted too. Since at issue is only the correlation among the residuals (leaving any correlations among the regressors alone), permuting the residuals (or the observations) will not generate the appropriate reference set for assessing the magnitude of the observed value of Moran's I.
  • In Cliff and Ord's original formulation of assumption N, new response variables are selected as random and independent drawings from a normal distribution. They argue this is the only workable assumption for regression residuals.
• For ordinary normal-based regression models, we have two options on how to proceed.
  1. We could for each observation generate a new observation by sampling from a normal distribution with mean , the value predicted by the fitted model, and variance , also estimated from the fitted model. Using the new set of response variables we then refit the model, obtain a new set of residuals, and calculate Moran's I for this new set of residuals.
  2. For normal-based models this process can be further simplified. Rather than obtaining new residuals by regressing a new set of (with a mean that varies in the fashion described by the model) on the old set of predictors, we can obtain a set of standard normal deviates and regress them on the predictors to obtain a new set of residuals.
• For a negative binomial regression model, option 2 is not available—the residuals do not have the same relationship with the response as they do in ordinary linear regression. So we're left with option 1. Option 1 applied to negative binomial regression uses the predicted mean value for each observation (obtained by applying the fitted function of R to the model) and the estimated size parameter as the parameters of a negative binomial distribution. The negative binomial random number generator of R, rnbinom, is then used to generate a value from this distribution. This is then done for each observation and a new regression model fit.
• The function below that I call residual.NB.MoranI (also found at http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/Moran residual.doc) carries out this protocol. It takes six arguments.
  1. The name of the data frame containing the data used in fitting the model.
  2. The name of the negative binomial regression object that contains the results of the regression fit using the actual data.
  3. The number of randomizations desired. Since the actual data values count as one randomization, this should be one less than the actual desired number.
  4. The name of the weight matrix containing the neighborhood information.
  5. The type of residual desired, "deviance" or "pearson", in quotes.
  6. An integer as a seed to initialize the random number stream. If a seed is not specified a randomly generated positive integer is used instead. The importance of specifying the seed yourself will be explained further below.

The function residual.NB.MoranI

# This function carries out randomization test of Moran I for negative binomial
# regression residuals. It uses glm.nb from the MASS package
# The five arguments are: argument 1 = dataset for fitting model
#   argument 2 = name of model output from NB regression for actual data
#   argument 3 = number of simulations desired, such as 999 for 1000 observations
#   argument 4 = the name of the weight matrix containing neighborhood information
#   argument 5 = the type of residual, "deviance" or "pearson", in quotes
#   argument 6 = numeric value for seed in generating simulated data

residual.NB.MoranI<-function(dataset,model,numsims,wtmat,restype,   
     seed=sample(1:1000000,1))
{
set.seed(seed)
#define moranI function
moranI.calc<-function(simres,wtres)
{
if (length(simres)>0)
{
simstd<- simres -mean(simres)
moranI<-t(simstd)%*% wtres %*% simstd/(t(simstd)%*% simstd)
} else moranI<-NA
moranI
}
#create matrix to hold results
moran.results<-NULL
flag<-0
i<-0
j<-0
#obtain model components
model$call->modelformula
RHS<-strsplit(as.character(modelformula[2]),'~')[[1]][2]
#this version only adds residuals if model converges
while(flag==0)
{
#generate new responses based on NB regression model
temp.y<-sapply(fitted(model),function(x) rnbinom(1,mu=x,size=model$theta))
#regress simulated data
glm.nb(as.formula(paste('temp.y', '~',RHS)),
   data=dataset,na.action=na.omit,control=glm.control(maxit=100))->simout1
#check if model converged. If no, do not add residuals.
if(simout1$converged) {
#obtain residuals and calculate Moran I
simres<-residuals(simout1,type=restype)
cur.moran<-moranI.calc(simres,wtmat)
i<-i+1
moran.results<-c(moran.results,cur.moran)
}
#next simulation
j<-j+1
#check if number of simulations requested is met
if(i==numsims) flag<-1
}
#finally add the actual values as the last row
res.actual<-residuals(model,type=restype)
last.moran<-moranI.calc(res.actual,wtmat)
moran.results<-c(moran.results,last.moran)
names(moran.results)<-c(as.character(1:numsims),"Actual")
moran.results
}

• The next function, plot.NB.MoranI uses the results of the residual.NB.MoranI function to plot the randomization distribution and calculate the randomization-based p-value. It takes three arguments.
  1. The name of the vector containing the output from the residual.NB.MoranI function.
  2. The name of the residual type, in quotes, only for purposes of labeling the graph.
  3. A 1 or a 2 to indicate a one-tailed or two-tailed p-value is desired.
• The p-value is calculated in the usual fashion for a randomization test. If the observed value of Moran's I is greater than the median of the randomization distribution, the program counts the number of randomization-based values of Moran's I that are greater than or equal to it (including the observed value itself) and divides by the total number of randomizations requested (adding 1 to include the observed value). If the observed value is below the median then randomization-based values that are smaller than or equal to the observed value are counted instead. If a two-tailed test is requested, the one-tailed number is simply doubled. If the one-tailed result is counting the wrong tail, e.g., you want an upper tail test and it reports a lower tail test, subtract the reported p-value from 1.

The function plot.NB.MoranI

#This function calculates p-value and plots randomization distribution
#arguments: argument 1 is the vector of simulated Moran I results
#   argument 2 is the name (in quotes) of the residual type for labeling graph
#   argument 3 should be 1 or 2 indicating a one- or two-tailed test
plot.NB.MoranI<-function(out.sims,restype,numtails=2)
{
#calculate p-value
cursim<-out.sims[!is.na(out.sims)]
actual.val<-cursim[length(cursim)]
#which tail is it in
medval<-median(cursim)
num.bigger<-ifelse(actual.val>medval,sum(cursim>=actual.val),
   sum(cursim<=actual.val))
#double one-tailed p-value for two-tailed value
p.val<-if(numtails==1) num.bigger/(length(cursim)) else
    2*num.bigger/(length(cursim))
morp<-out.sims[1:(length(out.sims)-1)]
zz<-density(morp)
plot.density(zz,main="Moran's I Residual Permutation Test",
    xlab="Reference Distribution",
    lwd=2,col=2,xlim=c(min(out.sims),max(out.sims)))
hist(morp,freq=F,add=T)
points(actual.val,0,pch=17,col=4,cex=1.5)
arrows(actual.val,0.4, actual.val,3,lwd=2,length=.09,code=1,col=4,lty=1)
tailtext<-if(numtails==1)   paste(restype, " residuals",
   ", # of simulations = ",
   prettyNum(length(cursim), big.mark=","),
   ", p = ",format.pval(p.val),' (one-tailed)',sep='') else 
   paste(restype, " residuals",", # of simulations = ",
   prettyNum(length(cursim), big.mark=","),
   ", p = ",format.pval(p.val),' (two-tailed)',sep='')
mtext(tailtext,line=0.5,side=3,cex=.8)
pretty.actual<-round(actual.val,3)
par(xpd=TRUE)
text(actual.val,3,"observed\nvalue ",col=4,cex=.8,pos=3)
par(xpd=FALSE)
}

#SAMPLE RUNS
#load MASS package
#library(MASS)
#carry out randomizations--999 translates into 1000 total observations
#test.sims<-residual.NB.MoranI(year5.disease,output.nb,999,wtmat5,"pearson",100)
# Calculate p-value and plot results. I do a one-tailed test.
#plot.NB.MoranI(test.sims,"pearson",numtails=1)

• Some implementation details.
  • Because of the way the new response variables are generated, it is possible to obtain simulated data set for which the negative binomial regression function, glm.nb, fails to converge. The program checks for this and if convergence did not occur, the solution is rejected and is not counted as one of the requested randomization results. As a result the program may run many more times than the requested number of randomizations.
  • Occasionally a particularly deadly combination of data values may be generated. I have found that the glm.nb function is prone to unceremoniously abort when faced with recalcitrant data. This also causes my function to stop running without saving the results from any of the models that had been successfully fit up to that point. If this happens there are two options.
  1. Restart using a different seed. This is less than ideal because a new seed may do no better.
  2. Use the same seed but reduce the number of randomizations in hopes that the last requested randomization occurs before the abortive run occurs. Having obtained at least some results with this seed, a new seed can be a tried to obtain the remaining number of desired randomizations. This step may need to be repeated more than once. The results can be then concatenated together to obtain the full number of randomizations.
  • Important: If option 2 is used and multiple runs are carried out it is important to drop the last observation of all runs but the last one. The residual.NB.MoranI program automatically appends the observed value of Moran's I to the simulated values as the last observation. If multiple runs are just concatenated into a single vector, then multiple copies of the observed value of Moran's I will also be included thus inflating the calculated p-value.
  • As an example, suppose three runs are required in which were requested 300, 200, and 499 randomizations. The actual length of the output vectors will be 301, 201, and 500. Only 1:300 of the observations from the first run, 1:200 observations from the second run, but all of the observations from the last run should be used.
• I demonstrate the program next. I start with a random seed of 10.

> residual.NB.MoranI(corals, model3, 999,wtmat5, "pearson", 10)->out1
Warning messages:
1: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
2: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
3: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
4: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
5: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
6: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
7: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
8: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,

• The warning messages are not a problem. They indicate that some of the runs failed to converge. These runs were discarded and do not affect the results. If we check the output we see it has the correct length.

> length(out1)
[1] 1000

• Finally I obtain the p-value and the randomization distribution (Fig. 6).

> plot.NB.MoranI(out1, "pearson", 1)

• Based on the graph we would not deem the observed value of Moran's I to be statistically significant at α = .05. If we were to repeat this for different seeds we would find p-values fluctuating in the range of .04 to .07.
• To understand what's going on we can plot the regression residuals versus geographic sector. Box plots are a standard way to display multiple distributions on the same plot. Since the sample size is small it's probably a good idea to show the raw data also. (For a discussion of this issue from an ecological perspective see Magnusson 2000.) I superimpose the jittered values of the observations (jittered to minimize overlap) on top of the box plots (Fig. 7).

> boxplot(residuals(model3, type='pearson') ~corals$SECTOR, outline=FALSE)
> points(jitter(as.numeric( corals$SECTOR)), residuals(model3, type='pearson'), pch=16, col=2 , cex=.8)
> mtext('Pearson residuals by sector', side=3, line=.5, cex=.9)

• The plot indicates why the spatial correlation is close to being statistically significant. Observe that the three largest residuals all occur in a single sector CA. Also, all but one of the seven residuals in each of sectors TO and SW are negative. Other than that there appears to be nothing unusual going on. Sectors CL and WH look interchangeable, and CB doesn't look too deviant given its small sample size.
• The analysis thus far has used Pearson residuals. Deviance residuals are generally preferred. I redo the analysis this time using deviance residuals. I use the same seed that was used with the Pearson residuals.

> residual.NB.MoranI(corals, model3, 999, wtmat5, "deviance", 10)->out2
Warning messages:
1: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
2: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
3: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
4: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
5: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
6: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
7: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
8: algorithm did not converge in: glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,
> length(out2)
[1] 1000
> plot.NB.MoranI(out2, "deviance", 1)
> boxplot(residuals(model3, type='deviance')~corals$SECTOR,outline=FALSE)
> points(jitter(as.numeric(corals$SECTOR)), residuals(model3, type='deviance'), pch=16, col=2 , cex=.8)
> mtext('Deviance residuals by sector', side=3, line=.5, cex=.9)

• The deviance residuals are better behaved. Based on these the evidence for spatial autocorrelation is weaker. If we plot the deviance residuals against sector we see that they are better behaved although the residuals from sectors TO and SW are still predominantly negative.

Fig. 8  Results for deviance residuals

Troubleshooting the randomization program

• In this section I discuss some of the problems alluded to above that can arise when constructing the randomization distribution. It's possible when generating new response values for the randomization distribution to obtain data that yield a model that doesn't converge. For example, suppose we elected to assess spatial correlation in the residuals of model2c fit above. Proceeding as before, we obtain the following.

> residual.NB.MoranI(corals,model2c,999,wtmat5,"pearson",10)->out1a
Error: NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message:
step size truncated due to divergence
> out1a
Error: Object "out1a" not found

• The "Error" message indicates a fatal error. If we check the contents of out1a it is empty. I cut the number of randomizations in half and try again.

> residual.NB.MoranI(corals,model2c,500,wtmat5,"pearson",10)->out1a
Error: NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message:
step size truncated due to divergence
> out1a
Error: Object "out1a" not found

• This one also aborted. I drop the request by 100. This also fails. When I drop the number of randomizations to 300 it works.

> residual.NB.MoranI(corals, model2c, 400, wtmat5,"pearson", 10)->out1a
Error: NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message:
step size truncated due to divergence
> residual.NB.MoranI(corals, model2c, 300, wtmat5, "pearson", 10)->out1a
> length(out1a)
[1] 301

• I next change the seed to 100 and try for the remaining randomizations.

> residual.NB.MoranI(corals, model2c, 699, wtmat5, "pearson", 100)->out1b
Warning messages:
1: iteration limit reached in: theta.ml(Y, mu, n, limit = control$maxit, trace = control$trace >
2: NaNs produced in: sqrt(1/i)
3: iteration limit reached in: theta.ml(Y, mu, n, limit = control$maxit, trace = control$trace >
4: NaNs produced in: sqrt(1/i)
> length(out1b)
[1] 700

• The "Warning messages:" look ominous, but actually is in reference to the iterations that failed to converge. These were rejected by the program and as a result were not added to the output vector. The output vector has 700 observations as desired.
• I next append the two sets of results being sure to drop the last observation from the first vector so that the observed value of Moran's I does not get counted twice.

> out<-c(out1a[1:300], out1b)

• Finally I can obtain the p-value and the randomization distribution as before (output not shown).

> plot.NB.MoranI(out, "pearson", 1)

Temporal autocorrelation (not done in class)

• As was demonstrated in lecture 4, reefs that are geographically close were also all sampled very close in time. Thus it would not be terribly surprisingly if the original response values exhibited a strong temporal correlation. Temporal correlation among regression residuals is a violation of stochastic indepedence and can affect the conclusions drawn from statistical inference. It's possible though for the temporal correlation to be "explained" by regression model. Thus the temporal correlation in prevalence could be accounted for by its relationship with WSSTA and coral cover.
• The presence or absence of significant temporal correlation in the regression residuals can be demonstrated by constructing autocorrelation plots in which correlation is plotted for different lag times. I begin by sorting the observations in the coral data frame by order of their sample date.

> library(date) #need the as.date function from this package
> ordered.disease<-corals[order(as.date(as.character(corals$DATE))),]

• Examining the data we see that the timing of the samples is fairly irregular.

> as.character(ordered.disease$DATE)
 [1] "7/24/2002"  "7/26/2002"  "7/28/2002"  "7/28/2002"  "7/30/2002"  "7/31/2002"
 [7] "7/31/2002"  "8/1/2002"   "8/2/2002"   "9/5/2002"   "9/6/2002"   "9/7/2002"
[13] "9/8/2002"   "9/9/2002"   "9/10/2002"  "9/11/2002"  "9/16/2002"  "9/17/2002"
[19] "9/18/2002"  "9/19/2002"  "10/24/2002" "10/25/2002" "10/26/2002" "10/28/2002"
[25] "11/1/2002"  "11/2/2002"  "11/3/2002"  "11/5/2002"  "11/7/2002"  "11/8/2002"
[31] "11/25/2002" "11/26/2002" "11/27/2002" "11/28/2002" "11/30/2002" "12/1/2002"
[37] "12/3/2002"  "12/4/2002"  "12/5/2002"  "3/7/2003"   "3/8/2003"   "3/10/2003"
[43] "3/12/2003"  "3/15/2003"  "3/18/2003"  "3/19/2003"  "3/21/2003"

• I next step through these observations calculating how far apart in time the observations are. I do this for varying apparent lags (how far apart observations are in the ordered list) each time calculating the true lag (the actual difference in date values).

> sapply(1:10,function(i) as.date(as.character(ordered.disease$DATE))[(1+i): (dim(ordered.disease)[1])]-
as.date(as.character(ordered.disease$DATE))[1: (dim(ordered.disease)[1]-i)]) ->the.lags

• The result is a list of true time lags. For example, the first component in the list created above lists the actual time separation between adjacent observations when ordered by date.

> the.lags[[1]]
[1]  2 2 0 2 1 0 1 1 34 1 1 1 1 1 1 5 1 1 1 35 1 1 2 4 1 1 2 2 1 17 1 1 1 2 1 2 1 1 92
[40] 1 2 2 3 3 1 2

• Observe that some observations were made on the same day but that most observations were made one day apart. The next list element in this variable computes the true date differences between paris of observations with apparent lag = 2 (those that are separated from each other by one other observation).

> the.lags[[2]]
[1]  4 2 2 3 1 1 2 35 35 2 2 2 2 2 6 6 2 2 36 36 2 3 6 5 2 3 4 3 18 18 2 2 3 3 3 3 2 93 93
[40] 3 4 5 6 4 3

• Notice there are still some observations in this vector whose true lag is one day. If we look at the last component of the list we can see what lags have been accounted for.

> min(the.lags[[10]])
[1] 14

• So, we can be sure we've located all pairs of observations whose true sample dates are within 13 days of each other. I tabulate all the true lags.

> table(unlist(the.lags))

 0  1  2  3  4  5  6  7  8   9  10  11  12  13  14  15  17  18  19  20  21  22  23  24  25
 2 27 28 22 18 16 14 14 14  11   9   9   6   5   4   1   1   2   2   3   2   3   5   4   5
26 27 28 29 30 31 32 34 35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
 5  4  5  4  3  3  3  1  3   6   8   8  10   9   8   7   5   5   7   8  10   6   3   2   1
52 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 119
 2  1  2  2  2  2  4  2  3   4   4   5   3   5   4   2   3   3   1

• So, two pairs of obseervations were made on the same day, 27 paris of observations were made one day apart, etc. Based on the reported numbers for lags of 8 or less we should have enough data to estimate the correlation coefficient. The function get.acf shown below and available at http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab7/ACF function.txt carries out the calculations. It locates all pairs of observations that have the same true lag and then uses them in calculating the empirical autocorrelation coefficient. It starts by standardizing the input variable as follows.

Using the standardized version it then calculates the autocorrelation at lag l, , as shown below (Pinheiro & Bates 2000, p. 227).

Here N(l) is the number of pairs that go into calculating the correlation while N(0) is the total number of observations that are included in the numerator at some point. M refers to the various different apparent lag groupings contained in the list variable the.lags.

The function get.acf

#This function calculates the autocorrelation function for a variable
# Arguments: argument 1 is the maximum lag to use in calculating the ACF
# argument 2 is the variable to use in calculating ACF
# argument 3 is the list variable containing the true lag values
get.acf<-function(maxlag,oldei,the.lags) {
 #standardize variable
 ei<-(oldei-mean(oldei))/sd(oldei)
 #initialize variables
 N.list<-length(ei) 
 ACF.list<-1
 total.sum<-0
 total.count<-0
 #cycle through all desired lags
 for (i in 1:maxlag) {
   all.guys<-NULL
   full.list<-NULL
  #cycle through lag lists looking for desired lag
    for (j in 1:length(the.lags)) {
       cur.lag<-the.lags[[j]]
   #calculate numerator terms of Morans I
       first.ei<-ei[((1+j):(length(the.lags[[j]])+j))[the.lags[[j]]==i]]
       second.ei<-ei[(1:(length(the.lags[[j]])))[the.lags[[j]]==i]]
       cur.term<-first.ei*second.ei
       full.list<-c(full.list,cur.term)
    #determine which observations were used for this lag
       cur.guys<-unique(c(((1+j):(length(the.lags[[j]])+j))[the.lags[[j]]==i],
       (1:(length(the.lags[[j]])))[the.lags[[j]]==i]))
       all.guys<-unique(c(all.guys,cur.guys))
       }
    #average sum of squares of all observations used
    bottom<-sum(ei[all.guys]^2)/length(ei[all.guys])
  #average lagged sum of squares
    top<-sum(full.list)/length(full.list)
    cur.ACF<-top/bottom
    cur.N<-length(full.list)
    ACF.list<-c(ACF.list,cur.ACF)
    N.list<-c(N.list,cur.N)
    }
 out.mat<-cbind(0:maxlag,ACF.list,N.list)
 colnames(out.mat)<-c('lag','ACF','N')
 out.mat
}

Fig. 9  Autocorrelation of the response

• I begin by trying out the function on the response variable, PREV_1.

> out.acf<-get.acf(10, ordered.disease$PREV_1, the.lags)
> out.acf

     lag         ACF  N
[1,]   0  1.00000000 47
[2,]   1  0.98290774 27
[3,]   2  0.60782175 28
[4,]   3  0.54641949 22
[5,]   4  0.35289774 18
[6,]   5  0.48385855 16
[7,]   6  0.24276398 14
[8,]   7  0.11626542 14
[9,]   8 -0.00919083 14
[10,]  9 -0.22915190 11
[11,] 10 -0.19324603  9

• Critical bounds for the autocorrelation are given by where is the standard normal quantile (Pinheiro & Bates 2000, p. 241). I plot the empirical autocorrelations along with the two-sided critical bounds (Fig. 9).

> plot(out.acf[,1], out.acf[,2], type='h', xlab='Lag', ylab='Autocorrelation', ylim=c(-1,1), col=4)
> abline(h=0,lty=1)
> lines(out.acf[,1], qnorm(.975)/sqrt(out.acf[,3]), lty=2, col=2)
> lines(out.acf[,1], -qnorm(.975)/sqrt(out.acf[,3]), lty=2, col=2)
> mtext('Autocorrelation of response', line=.5, side=3, font=2, cex=.9)

• The plot reveals that the empirical autocorrelations at the first three lags are significantly different from zero. I next turn to the residuals from the regression model.

> model3<-glm.nb(PREV_1 ~ WSSTA + CORAL_COVE + I(WSSTA^2), data = corals, control=glm.control(maxit=100))
> out.acf<-get.acf(10, residuals(model3, type='deviance'), the.lags)
> plot(out.acf[,1], out.acf[,2], type='h', xlab='Lag',ylab='Autocorrelation', ylim=c(-1,1), col=4)
> abline(h=0, lty=1)
> lines(out.acf[,1], qnorm(.975)/sqrt(out.acf[,3]), lty=2, col=2)
> lines(out.acf[,1], -qnorm(.975)/sqrt(out.acf[,3]), lty=2, col=2)
> mtext('Autocorrelation of residuals', line=.5, side=3, font=2, cex=.9)

• As we can see there is no evidence of any residual autocorrelation (Fig. 10). The original temporal correlation in the response variable has been entirely eliminated by the regression model.

Cited References

• Cliff, A. D. and J. K. Ord. 1981. Spatial Processes: Models & Applications. Pion Limited: London.
• Dunn, K. P. and G. K. Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics 5: 1–10.
• Magnusson, W. E. 2000. Error bars: are they the king's clothes? Bulletin of the Ecological Society of America 81: 147–150.
• Pinheiro, J. C. and D. M. Bates. 2000. Mixed-Effects Models in S and S-Plus. Springer-Verlag: New York.
• Venables, W. N. and B. D. Ripley. 2002. Modern Applied Statistics with S. Springer-Verlag: New York.

Course Home Page