Moth ExampleGeneralized Linear Models Case StudiesBret LargetDepartments of Botany and of StatisticsUniversity of Wisconsin—MadisonFebruary 26, 20081 / 13Poisson Regression Moth ExampleResearchers studied the effect of habitat on density of two differentmoth species.Researchers looked for moths along transects.Each transect was partitioned into different sections depending on thehabitat type.There are 8 different types of habitat.Counts of two moth species (labeled A and P) were obtained for eachpart of each transect.Under a Poisson model, we expect the counts to be proportional tothe length of the transect within each habitat.Moth Example 2 / 13Dotplot of DataMoth Density (per 10m)AP0 5 10 15 20●●Bank●●●● ● ●●● ● ●●● ●●Disturbed0 5 10 15 20●● ●● ●● ● ●●●●●●●●● ●●Lowerside●●●●●●●●●●●●NEsoakAP● ● ●●● ●NWsoak0 5 10 15 20● ●● ● ●● ●●●●●●●●SEsoak● ●●● ●●SWsoak0 5 10 15 20●●●●●●●●●●UppersideMoth Example 3 / 13Scatterplot of DataP Moth Density (per 10m)A Moth Density (per 10m)051015200 5 10 15●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●BankDisturbedLowersideNEsoakNWsoakSEsoakSWsoakUpperside●●●●●●●●Moth Example 4 / 13Model with an OffsetThe count data for species A comes from transects of differentlengths.Poisson process assumptions suggest that the expected count for asingle transect should be proportional to the length of the transect.We need to modify our typical Poisson regression model toaccomodate an exposure variable.Moth Example 5 / 13Model DerivationFor transect i:Iuiis the length in meters;Ihiis the habitat;Iθi= θ(hi) is the rate;Iyib e the count of species A;The expected count for the i th habitat is then:E[yi] = µi= uiθiMoth Example 6 / 13Derivation (cont.)We model θi= exp(Xiβ) where Xiincludes all predictors except forlength.In this case:Iθi= exp(β1) if hiis the first habitat;Iθi= exp(β1+ βk) if hiis the kth habitat, k > 1.We need E[yi] = µi= uiθi.With the exponential inverse link function, this is:µi= exp(log ui+ log(θi))=exp(log ui+ β1) if k = 1exp(log ui+ β1+ βk) if k > 1This is like regular Poisson regression with log u as a predictor, exceptwe want to force the estimated coefficient to be equal to one.Moth Example 7 / 13Fitting the Model> fitA = with(moths, glm(A ~ habitat, data = moths, family = poisson,+ offset = log(length)))> fitAq = with(moths, glm(A ~ habitat, data = moths, family = quasipoisson,+ offset = log(length)))> fitP = with(moths, glm(P ~ habitat, data = moths, family = poisson,+ offset = log(length)))> fitPq = with(moths, glm(P ~ habitat, data = moths, family = quasipoisson,+ offset = log(length)))Moth Example 8 / 13Poisson Regression for A> display(fitA)glm(formula = A ~ habitat, family = poisson, data = moths, offset = log(length))coef.est coef.se(Intercept) -18.35 1275.75habitatDisturbed 16.53 1275.75habitatLowerside 16.81 1275.75habitatNEsoak 15.45 1275.75habitatNWsoak 18.44 1275.75habitatSEsoak 16.70 1275.75habitatSWsoak 16.59 1275.75habitatUpperside 14.82 1275.75---n = 41, k = 8residual deviance = 180.3, null deviance = 499.0 (difference = 318.7)Moth Example 9 / 13QuasiPoisson Regression for A> display(fitAq)glm(formula = A ~ habitat, family = quasipoisson, data = moths,offset = log(length))coef.est coef.se(Intercept) -18.35 3793.81habitatDisturbed 16.53 3793.81habitatLowerside 16.81 3793.81habitatNEsoak 15.45 3793.81habitatNWsoak 18.44 3793.81habitatSEsoak 16.70 3793.81habitatSWsoak 16.59 3793.81habitatUpperside 14.82 3793.81---n = 41, k = 8residual deviance = 180.3, null deviance = 499.0 (difference = 318.7)overdispersion parameter = 8.8Moth Example 10 / 13Poisson Regression for P> display(fitP)glm(formula = P ~ habitat, family = poisson, data = moths, offset = log(length))coef.est coef.se(Intercept) -1.66 0.50habitatDisturbed 1.26 0.53habitatLowerside -0.76 0.56habitatNEsoak -1.24 0.57habitatNWsoak 0.43 0.55habitatSEsoak -1.81 0.65habitatSWsoak 0.78 0.52habitatUpperside -3.12 0.61---n = 41, k = 8residual deviance = 92.6, null deviance = 389.5 (difference = 296.9)Moth Example 11 / 13QuasiPoisson Regression for P> display(fitPq)glm(formula = P ~ habitat, family = quasipoisson, data = moths,offset = log(length))coef.est coef.se(Intercept) -1.66 1.07habitatDisturbed 1.26 1.13habitatLowerside -0.76 1.19habitatNEsoak -1.24 1.21habitatNWsoak 0.43 1.18habitatSEsoak -1.81 1.38habitatSWsoak 0.78 1.11habitatUpperside -3.12 1.31---n = 41, k = 8residual deviance = 92.6, null deviance = 389.5 (difference = 296.9)overdispersion parameter = 4.6Moth Example 12 / 13MoreMore live with R and on BoardMoth Example 13 /
View Full Document
Unlocking...