Search this blog

A function of Pearson/Deviance goodness-of-fit for (generalized) linear regression

After an analysis of (generalized) linear regression, a Pearson/Deviance goodness-of-fit test is useful to test if the model fitted reasonably. I wrote a simple function in R to do it.

Source Code

#    Copyright 2011 Chen-Pan Liao 
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    any later version.
#    
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#    
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

# Usage example:
#   model.name <- glm (...)
#   lm.fit.test(model.name)

lm.fit.test <- function (my.model) {
  df <- my.model$df.residual;
  chisq.pearson <- sum(resid(my.model , type="pearson")^2);
  chisq.deviance <- sum(resid(my.model , type="deviance")^2);
  p.pearson <- pchisq( chisq.pearson , df , lower.tail=F );
  p.deviance <- pchisq( chisq.deviance , df , lower.tail=F );
  ratio.pearson <- chisq.pearson / df;
  ratio.deviance <- chisq.deviance / df;
  cat(
    "Pearson chisq = " , chisq.pearson ,
    ", df = " , df ,
    ", p = " , p.pearson ,
    ", chisq / df = " , chisq.pearson/df ,
    ".\n" ,
    "Deviance chisq = " , chisq.deviance ,
    ", df = " , df ,
    ", p = " , p.deviance ,
    ", chisq / df = " , chisq.deviance/df ,
    ".\n"
  );
}

Usage and example

After defining a glm() model as a custom name, copy the source code into R interpreter and then input lm.fit.test(your.model.name) where your.model.name is the custom glm model name. See the fallowing example

y <- c(rpois(20,1) , rpois(20,2));
x <- gl(2 , 20);
my.model <- glm(y~x , family=poisson);
summary(my.model);
lm.fit.test(my.model);