Skip to content

Instantly share code, notes, and snippets.

# hinkelman/CountGLM.R Created Apr 15, 2013

 n = 100 # No. of samples per treatment mean.A = 10 # Mean count for treatment A mean.B = 5 # Mean count for treatment B nd = data.frame(Trt = c("A", "B")) # Used in predict( ) function # Generate example data based on Poisson distribution (mean = variance) data.pois = data.frame(Trt = c(rep("A", n), rep("B", n)), Response = c(rpois(n, mean.A), rpois(n, mean.B)) ) # Generalized linear model with Poisson distribution of errors model.pois = glm(Response ~ Trt, data = data.pois, family = poisson) summary(model.pois) # Test for goodness of fit of model; compare residual deviance with df (from summary output) # If test is significant, model fit is not good 1 - pchisq(summary(model.pois)\$deviance, summary(model.pois)\$df.residual ) # Predict mean counts and standard errors cbind(nd, Mean = predict(model.pois, newdata=nd, type="response"), SE = predict(model.pois, newdata=nd, type="response", se.fit=T)\$se.fit ) #-------------------------------------------------------------------------- # Generate example data based on negative binomial distribution (mean < variance) library(MASS) data.nb = data.frame(Trt = c(rep("A", n), rep("B", n)), Response=c(rnegbin(n, mean.A, 5), rnegbin(n, mean.B, 5)) ) # Generalized linear model with Poisson distribution of errors model.pois.2 = glm(Response ~ Trt, data = data.nb, family = poisson) summary(model.pois.2) # Test for goodness of fit of model; compare residual deviance with df (from summary output) 1 - pchisq(summary(model.pois.2)\$deviance, summary(model.pois.2)\$df.residual ) # Predict mean counts and standard errors cbind(nd, Mean = predict(model.pois.2, newdata = nd, type = "response"), SE = predict(model.pois.2, newdata = nd, type = "response", se.fit = T)\$se.fit ) # Generalized linear model with negative binomial distribution of errors model.nb = glm.nb(Response ~ Trt, data = data.nb) summary(model.nb) # Test for goodness of fit of model; compare residual deviance with df (from summary output) 1 - pchisq(summary(model.nb)\$deviance, summary(model.nb)\$df.residual ) # Predict mean counts and standard errors cbind(nd, Mean = predict(model.nb, newdata = nd, type = "response"), SE = predict(model.nb, newdata = nd, type="response", se.fit = T)\$se.fit ) #-------------------------------------------------------------------------- # Generate example data based on zero-inflated Poisson distribution library(VGAM) data.zip = data.frame(Trt = c(rep("A", n), rep("B", n)), Response = c(rzipois(n, mean.A, 0.2), rzipois(n, mean.B, 0.2)) ) # Generalized linear model with Poisson distribution of errors model.pois.3 = glm(Response ~ Trt, data = data.zip, family = poisson) summary(model.pois.3) # Test for goodness of fit of model; compare residual deviance with df (from summary output) 1 - pchisq(summary(model.pois.3)\$deviance, summary(model.pois.3)\$df.residual ) # Predict mean counts and standard errors cbind(nd, Mean = predict(model.pois.3, newdata = nd, type = "response"), SE = predict(model.pois.3, newdata = nd, type = "response", se.fit = T)\$se.fit ) # Generalized linear model with negative binomial distribution of errors model.nb.2 = glm.nb(Response ~ Trt, data = data.zip) summary(model.nb.2) # Test for goodness of fit of model; compare residual deviance with df (from summary output) 1 - pchisq(summary(model.nb.2)\$deviance, summary(model.nb.2)\$df.residual ) # Predict mean counts and standard errors cbind(nd, Mean = predict(model.nb.2, newdata = nd, type = "response"), SE = predict(model.nb.2, newdata = nd, type = "response", se.fit = T)\$se.fit ) # Zero-inflated Poisson model library(pscl) # Test for treatment effect on counts model.zip = zeroinfl(Response ~ Trt|1, data = data.zip) summary(model.zip) # Predict mean counts and zeros cbind(nd, Count = predict(model.zip, newdata = nd, type = "count"), Zero = predict(model.zip, newdata = nd, type = "zero") ) # Test for treatment effect on counts and zeros model.zip.2 = zeroinfl(Response ~ Trt|Trt, data = data.zip) summary(model.zip.2) # Predict mean counts and zeros cbind(nd, Count = predict(model.zip.2, newdata = nd, type = "count"), Zero = predict(model.zip.2, newdata = nd, type = "zero") ) # Zero-inflated negative binomial model model.zip.3 = zeroinfl(Response ~ Trt|1, data = data.zip, dist = "negbin") summary(model.zip.3)
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.