Unformatted text preview:

Stat 544, Lecture 14 1'&$%More AboutGLIM DiagnosticsLast time, we began to discuss diagnostics forgeneralized linear models (GLIMs). These diagnosticsare motivated by the fact that we fit GLIMs byiteratively reweighted least squares,β(t+1)=“XTWX”−1XTWz,whereW = Diag"Var(yi)„∂ηi∂μi«2#−1is the matrix of weights andz = η +„∂η∂μ«(y − μ)is the working variate.Thus far, we have discussed• plotting the residuals (either Pearson ordeviance) versus the fitted values of the linearStat 544, Lecture 14 2'&$%predictor η. This helps us to look for outliers,overdispersion or an incorrect variance function.• plotting absolute residuals versus fitted values.This reveals problems with the variance function.Now we’ll discuss how to check the link function andhow to look for influential observations.Checking the link function. The simplest way tocheck the link function is to plot the final value of theworking variate z against the linear predictor η. Theplot should resemble a straight line. Curvature in thisplot suggests that the link function is not appropriate.For binary data, this plot is uninformative!Hinkley (1985) suggests the following test: Afterfitting the GLIM, calculate η2and then try adding itto the model as a covariate. If this covariate issignificant, then there is a problem. Significance inthis test could be caused by a wrong link function, awrong scale for one or more predictors, or both.Leverage. In ordinary regression, the diagonalelements hiiof the “hat matrix”H = X(XTX)−1XTStat 544, Lecture 14 3'&$%are called the leverage values. Because H isidempotent,H2= H,its rank equals its trace, and thus the average value ofhiiis be p/N. Values greater than, say, 3p/N areconsidered to have high leverage and may have highinfluence. For a GLIM, the hat matrix is obtained byreplacing X with W1/2X,H = W1/2X(XTWX)−1XTW1/2.An observation whose x is far from the centroid mightnot have high leverage if its weight is small.Plotting the residuals versus the leverage values ishelpful for detecting outliers and potentiallyinfluential observations at the same time.Example. Let’s apply these diagnostic techniques toa data example from McCullagh and Nelder (1989).This table shows the proportion of carrots damaged inan insecticide experiment carried out in three blocks.Stat 544, Lecture 14 4'&$%Blocklog-dose 1 2 31.52 10/35 17/38 10/341.64 16/42 10/40 10/381.76 8/50 8/33 5/361.88 6/42 8/39 3/352.00 9/35 5/47 2/492.12 9/42 17/42 1/402.24 1/32 6/35 3/222.36 2/28 4/35 2/31Before applying a GLIM, let’s consider what we woulddo if the response was continuous and heteroscedastic.The classical ANOVA table would look like this:Source SS dfBlock SSB2Treatment SST7Error SSB×T14We’d obtain an error term by assuming that there areno interactions, i.e. by assuming a constant treatmenteffect across blocks. If we did not assume this, wecould not test the significance of block or treatment.This ANOVA is equivalent to fitting a linearregression model with an intercept, two dummies todistinguish among the blocks, and seven dummies toStat 544, Lecture 14 5'&$%distinguish among levels of the treatment. Becausethe treatment levels are numeric, it makes sense toinvestigate simpler models that assume a simplerelationship (e.g. linear or quadratic) between theresponse and log-dose.Now let’s analyze the data by a GLIM. It makes senseto start out by assuming a binomial response and alogit link. We ought to consider adding block effectsby defining dummy indicators to distinguish amongthe blocks. The most complex model for thetreatment would have dummy indicators todistinguish among the treatment levels. But we couldalso try simpler models that assume simple (e.g.linear or quadratic) effects.Plotting the empirical logits. Because theresponse is a proportion, we should use a binomialmodel. Lacking any prior information about how theproportion might vary as a function of log-dose, wemight want to try the logistic link first, because it isthe canonical link and the most common one for abinomial model. To see how the logit-probability isrelated to log-dose, we may simply plot the empiricallogits versus log-dose.Stat 544, Lecture 14 6'&$%ooooooooooooooooooooooooThis plot shows a decreasing, possibly linear trendwhich suggests that a logit model with linear term forlog-dose might work.Let’s fit the logistic model with two dummy variablesfor block and a linear effect for log-dose. The resultsfrom SAS are shown below.Deviance and Pearson Goodness-of-Fit StatisticsCriterion DF Value Value/DF Pr > ChiSqDeviance 20 39.9757 1.9988 0.0050Pearson 20 41.1746 2.0587 0.0035Type III Analysis of EffectsWaldEffect DF Chi-Square Pr > ChiSqblock 2 13.9181 0.0009logdose 1 27.9323 <.0001Stat 544, Lecture 14 7'&$%Analysis of Maximum Likelihood EstimatesStandard WaldParameter DF Estimate Error Chi-Square Pr > ChiSqIntercept 1 1.4802 0.6562 5.0884 0.0241block 1 1 0.5424 0.2318 5.4749 0.0193block 2 1 0.8432 0.2260 13.9175 0.0002logdose 1 -1.8174 0.3439 27.9323 <.0001The effects of block and log-dose are both highlysignificant, but the model doesn’t fit,P (χ220≥ 41.2) = 0.0035. The possible reasons for lackof fit are:• the effect of log-dose might not be linear (but wesaw no evidence of that in the empirical logitplot);• there may be interactions between block andlog-dose; and• overdispersion.We probably do not want to fit a model withinteractions. Recall that the classical ANOVA modelfor a dataset like this assumes no interactions.Overdispersion is a clear possibility.Checking for overdispersion. To check foroverdispersion, we should fit a maximal model. ByStat 544, Lecture 14 8'&$%analogy to the classical ANOVA, we could say that anatural choice for the maximal model is an additivemodel with an intercept, two dummies to distinguishamong the blocks, and seven dummies to distinguishamong the levels of treatment. This model has nointeractions, but it does not assume any particularfunctional form (e.g. linear) between the response andlog-dose. Here’s SAS code to fit this model and togenerate two diagnostic plots.options nocenter nodate nonumber linesize=72;data carrot;input logdose block y n;cards;1.52 1 10 351.52 2 17 381.52 3 10 341.64 1 16 42- -lines omitted--2.36 2 4 352.36 3 2 31;proc logist data=carrot;class block (order=data param=ref ref=last);class logdose (order=data param=ref ref=first);model y/n = block logdose / scale=none;output out=out1 xbeta=xb


View Full Document

PSU STAT 544 - GLIM Diagnostics

Download GLIM Diagnostics
Our administrator received your request to download this document. We will send you the file to your email shortly.
Loading Unlocking...
Login

Join to view GLIM Diagnostics and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view GLIM Diagnostics 2 2 and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?