Create a gist now

Instantly share code, notes, and snippets.

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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment