Generalized Linear Modelling With Mathematica

Department of Statistical and Actuarial Sciences
Student : Paula Lau
Supervisor : A.I.McLeod
1998 Winter Term

Introduction

Regression analysis is very commonly used in data analysis to study the relationship between a response variable and one or more explanatory variables. To use regression analysis, we have to make sure that the underlying assumptions are meet. These assumptions includes the data has to be continuous and normally distributed, observations has to be uncorrelated and homogeneous variance across levels of the independent variable.  Least square is used for parameter estimation in linear models.

However, we often encounter practical problems where the outcome variable is not continuous or that the Gaussian error assumption is not valid.  For example, in medical trials with subjects being patients who had heart attack, we want to know how different drug treatments affect the chances of patients suffering recurrent heart attack at the end of the trial period. In this case the responses are either Yes or No, so the conditional distribution of the outcome variable and the error both follow a binomial distribution. The assumptions for linear regression are clearly violated and ordinary least-squares does not apply.

In situations like this, we need to fit a Generalized Linear model, GLM allows fitting linear regression model with general family of error distribution, it can be fitted for the data that has probability function a member of the exponential family.  The normal distribution of errors becomes a special case of this family of linear models.  Members in the expoential family include binomial, poisson and Gamma, exponential, negative binomial and others.

Definition

Model

Generalized Linear models are an extension of classical linear models.
Given a set of observations (Subscript[y, i]) (i=1, 2,.., n)  and covariates (Subscript[x, i]) (i=1, 2,., p)
and E((Subscript[y, i]))=(Subscript[µ, i]) 
Generalized Linear Models  consist of 3 parts :

1. linear predictor : (Subscript[(η ), i])= \(TraditionalForm\`Σ\+\(i = 1\)\%p\) (Subscript[x, i]) (Subscript[β, i])  
2. E((Subscript[y, i])) =(Subscript[µ, i])((Subscript[β, 1]), (Subscript[β, 2]),.., (Subscript[β, p]))
3. link : (Subscript[(η ), i])= g((Subscript[µ, i]))

where g(.) is called the link function.  The link function could be any monotonic differentiable function that maps the responses interval on to the whole real line (-∞, +∞).  Canonical links are commonly used since it has desirable statistical properties. We will discuss how we get canonical in the next section.
The canonical links for the most commonly used distributions are :

Normal : (Subscript[(η ), i]) = (Subscript[µ, i])
Poisson : (Subscript[(η ), i]) = log (Subscript[µ, i])
Binomial : (Subscript[(η ), i]) = log (((Subscript[µ, i]))/(1-(Subscript[µ, i])))
Gamma : (Subscript[(η ), i]) = (1/(Subscript[µ, i]))

Therefore, classical linear model with Normal errors becomes a special case of GLM with an identity link.   In GLM, we still assume observations independent of each other, the assumption on normal error distribution is no longer required.

Transformation as an Alternative to GLM

Data transformation is sometimes used as an alternative to Generalized linear model  to stabilize the variance, to improve the symmetry of the density and to produce additivity of the systematic effects.  We can then apply ordinary least squares to the transformed data.

For example, the variance-stabilizing transformation to binary data is arcsin; square root transformation  is used for poisson data and log transformation is used for Gamma data.

The drawback with data transformation is that a transformation for stability of variance is often different from the transformation to approximate Normality of errors and that for additivity.  So it is difficult to choose the best transformation for the data in hand.  With the introduction of Generalized Linear Models, Normality and Constancy of variance are no longer required.  Thus scaling problems are greatly reduced.

Although both methods lead to similar conclusions, the Generalized Linear Model has found to give a better fit in some cases.  We have tried the transformation approach for car insurance data which is fitted Generalized linear model with gamma error by McCullagh & Nelder and found that the data support the reciprocal model better.

Maximum Likelihood

The parameter estimates of Generalized Linear models are obtained by maximizing the likelihood of the parameters of the data observed, the maximum likelihood estimates are the same as the least square estimates when the errors are normal.  The log likelihood is defined as

l(μ; y) = (Underscript[Σ, i]) log(Subscript[f, i])(y, θ)

where f(y,θ) is the density function for the observation y given the parameter θ

For observations Y has a distribution in the exponential family taking the form

(Subscript[f, y]) (y; θ, ϕ) = 
      exp {(((yθ - b (θ)))/(a (ϕ))) +  
            c (y, ϕ)}

where θ is known as the canonical parameter. Canonical links occurs when θ = η.  ϕ is called the dispersion parameter.
The mean and variance of Y can be derived easily, so we have

E (Y) = µ  = (b^′) (θ)var (Y)=(b^″)(θ) a(ϕ)

(b^″)(θ) is known as the variance function and will be written V(µ).

The likelihood function has the form

l (θ, y) = 
      (((yθ - b (θ)))/(a (ϕ)))  +  
        c (y, ϕ)

Exact Likelihood

The usual method of estimating Generalized Linear models parameters is approximate maximum likelihood.  Since explicit solution to the minimization doesn't exist, most statistical packages like Splus and SAS use an iterated reweighted least square algorithm to get the approximate maximum likelihood estimates.  
Although some exact algorithms were developed later on, they tend to be rather computing intensive. So they were never considered to be a practical option due to the limitation in computing power.  
Until recently, with easy access to inexpensive computing power, we can make exact maximum likelihood computationally feasible.  
With Mathematica, we take advantage of the build-in function FindMin[] to obtain the exact maximum likelihood estimates as suggested by Currie.

Deviance and Deviance Residuals

The Scaled Deviance is a linear function of the log likelihood  l(μ;y), namely

D (μ; y) = 2 l (y; y) - 2 l (μ; y)

This is a measure of the difference between the log likelihood of the fitted model and that of the saturated model.  Saturated model is the model with fitted values equal to the observed data, so it has maximum achievable log likelihood.

The forms of deviances of the most commonly used distributions are:

Normal	Σ(((y-\(µ\&^\)))^2)
Poisson	2Σ{\(y\ \(log(y\/µ\&^)\)\)-(y-\(µ\&^\))}
binomial	Σ{\(y\ \(log(y\/µ\&^)\)\)+(m - y)log( ((m - y)/(m-\(µ\&^\))))}
gamma		Σ(\(-\ \(log(y\/µ\&^)\)\)+((y-\(µ\&^\))/(\(µ\&^\))))

Each observation contributes a quantity (Subscript[d, i]) to the deviance, so that  Σ(Subscript[d, i]) = D . 
Hence, the deviance residual is defined as

\(\(	r\_D = sign \((y - µ)\) \@d\_i\)\)

and Σ((Subscript[r, D])^2) = D.
Deviance residuals are used for diagnostic checks. The goodness-of-fit test of Generalized Linear model is also based on deviance.  These will be discussed in more details in the Model Selection Section.

Diagnostic Checks

Diagnostic checks for Generalized linear models are based on deviance residuals that are approximately normally and independently distributed under the hypothesis that the model is correct. 
The following diagnostic plots are included in the output of the  ExactGLM function.

1. Residuals vs Fitted Values
We expect to see residuals lies randomly about the line y=0 and doesn't show any patterns, otherwise it may suggest a lack of fit.

2. Sqrt of abs(Residuals) vs Predictions
This is the spread level plot and it can be used to detect if the variance changes with the level of the data.  An increasing pattern suggests the existence of monotone spread.

3. Response vs Fitted Values
For a good fit, we expect to see the data lies closely to the line y=x and doesn't show any patterns.

4. Normal QQplot of Std. Residual
It is a plot of the residuals quantiles against the normal N(0,1) quantiles.  If the residuals lies along or close to a straight line, it is well approximated by the normally distribution.

Examples

Binomial(grouped)

Fabric failure data : Varying number of fabric type specimens were exposed to load(x) in lb/(in^2).  And the number of specimens that failed is recorded.

fabric = {{5, 600, 13}, {35, 500, 95}, {70, 600, 189}, {80, 300, 
          95}, {90, 300, 130}};

ExactGLM[fabric, family = "binomial", {1, load}, {load}, 
      plot = True]

Binomial(ungrouped)

In an observation study, subjects were exposed to varying levels of rate and volume of air inspired on a transient vaso restriction on the skin of the digits.  Only the occurrence or nonoccurrence of the restriction was observed.

\(\(volrate = {{3.7, 0.825, 1, 1}, {3.5, 1.09, 1, 1}, {1.25, 2.5, 1, 1}, {
          0.75, 1.5, 1, 1}, {0.8, 3.2, 1, 1}, {0.7, 3.5, 1, 1}, {0.6, 0.75, 
          1, 0}, {1.1, 1.7, 1, 0}, {0.9, 0.75, 1, 0}, 

        		{0.9, 0.45, 1, 0}, {0.8, 0.57, 1, 0}, {0.55, 2.75, 1, 0}, {0.6, 
          3, 1, 0}, {1.4, 2.33, 1, 1}, {0.75, 3.75, 1, 1}, {2.3, 1.64, 1, 
          1}, {3.2, 1.6, 1, 1}, {0.85, 1.415, 1, 1}, {1.7, 1.06, 1, 0}, {1.8, 
          1.8, 1, 1}, {0.4, 2, 1, 0}, {0.95, 1.36, 1, 0}, {1.35, 1.35, 1, 
          0}, {1.5, 1.36, 1, 0}, {1.6, 1.78, 1, 1}, {0.6, 1.5, 1, 0}, {1.8, 
          1.5, 1, 1}, {0.95, 1.9, 1, 0}, 

        		{1.9, 0.95, 1, 1}, {1.60, 0.4, 1, 0}, {2.7, 0.75, 1, 1}, {2.35, 
          0.03, 1, 0}, {1.1, 1.83, 1, 0}, {1.1, 2.2, 1, 1}, {1.2, 2.0, 1, 
          1}, {0.8, 3.33, 1, 1}, {0.95, 1.9, 1, 0}, {0.75, 1.9, 1, 0}, {1.3, 
          1.625, 1, 1}
	}; \)\)

ExactGLM[volrate, 
      family = "binomial", {1, volume, rate}, {volume, rate}, plot = True]

Poisson

data2 = {{2, 1, 6}, {3, 2, 1}, {3, 3, 4}, {3, 3, 2}, {4, 2, 6}, {3, 
          1, 8}, {2, 1, 4}};

ExactGLM[data2, family = "poisson", {1, x, y}, {x, y}, plot = True]

Gamma

data1 = {{"b", 6}, {"c", 1}, {"c", 4}, {"c", 2}, {
          "d", 6}, {"c", 8}, {"b", 4}};

ExactGLM[data1, family = "gamma", {1, cat}, {cat}, plot = True]

Normal

data1 = {{"b", 6}, {"c", 1}, {"c", 4}, {"c", 2}, {
          "d", 6}, {"c", 8}, {"b", 4}};

ExactGLM[data1, family = "normal", {1, cat}, {cat}, plot = True]

Model Selection

Model Selection is the process of selecting the variables that result in the most parsimonious model that best explains the data.  Deviance is used in model selection process in Generalized Linear models like the way residual sum of squares are used in classical linear models.

To assess the significance of an independent variable, we need to balance between the goals of parsimony and improvement of fit.  We can apply the likelihood ratio test to resolve this problem.  Likelihood ratio test refers to comparing the value of deviance (D) with and without the variable in the model.

G=D(model without the variable)-D(model with the variable)

Under (Subscript[ℋ, 0]) : coefficient of variable is zero, G follow a chi-square distribution with 1 degree of freedom.  Rejection of  Null Hypothesis means that the variable is significant and  inclusion of the variable can improve the model.

The same approach can be used to test the significance of any subsets of variables to help us decide whether to choose a full model or a reduced model. So for a subset (Subscript[β, a]) with k<p variables, where (Subscript[β, c])  =  [(Subscript[β, a]) , (Subscript[β, b])] is the set of coefficients in our full model.  Under (Subscript[ℋ, 0]) : (Subscript[β, a])=0.  The test statistic is G=D(model with (Subscript[β, a])=0)-D(full model), G∼(𝒳^2)(k).   Rejection of Null Hypothesis means that at least one of the k coefficients is different from zero.

To test the overall significance of a model, we can compare the full model with the null model using G=D(null model)-D(full model) where G ∼(𝒳^2)(p).

After fitting the model, we would like to access the goodness-of-fit of the model, in other words, we want to know how effective the model is in describing the data.  Goodness-of-fit of GLM is measured using the Deviance as the test statistic.  Under the assumption that the fitted model is correct, D is supposed to be (𝒳^2)(n-(p+1)).

Example : To demonstrate model selection, we compare two alternative models with the full model we fitted.
volmodel : model with volume as the only explanatory variable

volmodel = 
      ExactGLM[volrate, 
        family = "binomial", {1, volume}, {volume,  rate}, plot = False]

ratemodel : model with rate as the only explanatory variable

ratemodel = 
      ExactGLM[volrate, family = "binomial", {1, rate}, {volume, rate}, 
        plot = False]

The deviance for the three models are computed :

\(voldev = 46.9894; 
ratedev = 49.6548\ ; 
fulldev = 29.7724; \)

To test the null hypothesis that rate is insignificant in the model, we compute the significance level of G=voldev-fulldev

1 - CDF[ChiSquareDistribution[1], voldev - fulldev]

Since p-value<0.05, there is strong evidence against the null hypothesis, thus rate should be included in the model.

To Test the null hypothesis that volume is insignificant in the model, we compute the significance level of G=voldev-fulldev

1 - CDF[ChiSquareDistribution[1], ratedev - fulldev]

Since p-value<0.05, there is strong evidence against the null hypothesis, thus volume should be included in the model.  So we conclude that both rate and volume are significant and should be included in the model.

The significance level of (𝒳^2) Goodness-of-fit test is included in the standard output of the ExactGLM function, the p-value of the full model is 0.758427, so there is no evidence against the statistical adequacy of the full model.

Program

Mathematica is becoming increasingly popular among statisticians due to its ability to perform symbolic computations.  However, the range of built-in function to perform statistical data analysis is still limited comparing with other statistical software.

The function  ExactGLM is written that fits  the Generalized Linear models, the algorithm takes advantage of the build-in FindMin[] function to obtain the exact parameter estimates.

References

Iain D. Currie (1994) Maximum Likelihood Estimation and Mathematica. Appl. Statist. (1995) 44, No. 3, pp. 379-394,
P. McCullagh and J.A. Nelder FRS (1989) Generalized Linear Models, 2nd Edition. Chapman and Hall
David W. Hosmer Stanley Lemeshow (1989) Applied Logistic Regression. John Wiley & Sons
R.H. Myers (2nd Edition) Classical and Modern Regression with Applications. PWS publishers.

Implementation of Exact Maximum Likehood Generalized Linear Model

η =  \(TraditionalForm\`Σ\+1\%p\) (Subscript[x, j])(Subscript[β, j])
              μ = g(η)
              g: canonical link function

set up the package context, included any imports

Get["Statistics`DiscreteDistributions`"]Get["Statistics`ContinuousDistributions`"]

usage messages for the exported functions and the context itself

ExactGLM :: usage = 
      "ExactGLM[data,family,fun,var] fits a Generalized linear model of \
family family to a list of data as a linear combination of the functions funs \
of variables vars. Canonical links are used as the link functions.";

set up private context

Begin["`Private`"]

MakeDummy

Create Dummy variables for  categorical independent variables, the function returns a Design Matrix with Categorical variables replaced by Dummy variables.

\(MakeDummy[Matrix_, funs_] := 
      Module[{k1, check, i, j, catvar, unique, \ level, varname, 
          newfuns = {}, NumericMatrix = {}}, 

        k1 = \(Dimensions[Matrix]\)[\([2]\)]; 

        check = Map[NumberQ, Matrix[\([1]\)]]; 

        Do[catvar = \(Transpose[Matrix]\)[\([i]\)]; 
			
          varname = funs[\([i]\)]; \ 
		\ \ 
          If[Not[check[\([i]\)]], 
				\ \ 
             (*\ do\ the\ following\ for\ categorical\ variable\ *) 

            				\ \ \ \ unique = Union[catvar, catvar]; 

            			\ \ \ \ \ \ \ level = Length[unique]; 
				\ \ \ \ 
            Do[
					\ \ \ \ \ \ \ \ 
              If[level == 2, \ 
						\ \ \ \ \ \ \ \ \ \ \ \ \ 
                newfuns = Append[newfuns, varname], 

                						\ \ \ \ \ \ \ \ \ \ \ \ \ 
                newfuns = \ 
                  Append[newfuns, 
                    Symbol[\ 
                      StringJoin[ToString[varname], 
                        ToString[unique[\([j + 1]\)]]]]]]; 

              					\ \ \ \ \ \ 
              NumericMatrix = 
                Append[NumericMatrix, 
                  catvar /. 
                    Join[{unique[\([j + 1]\)] -> 1}, \ \ 
                      Table[unique[\([i]\)] -> 0, {i, level}]]], \ \ 

              					\ {j, 1, level - 1}], 
				\ \ 
             (*\ do\ the\ following\ for\ non - categorical\ variable\ *) 

            				\ \ \ newfuns = Append[newfuns, varname]; 

            			\ \ \ \ \ \ NumericMatrix = Append[NumericMatrix, catvar]], 
          
{i, 1, k1}]; 
{Transpose[NumericMatrix], newfuns}
]\)

Link Functions

Canonical Links are used as the link functions

NORMAL

\(normal[y_, nullη_, η_, n_] := 
      Module[{nullμ, μ, k, nullprob, prob}, 
nullμ = nullη; 
        
μ = η; 

        nullprob = 
          Table[PDF[NormalDistribution[nullμ[\([k]\)], σ], 
              y[\([k]\)]], {k, 1, n}]; 

        prob = Table[
            PDF[NormalDistribution[μ[\([k]\)], σ], y[\([k]\)]], {
              k, 1, n}]; 
		{nullprob, prob}
	]\)

BINOMIAL

\(binomial[y_, nullη_, η_, m_, n_] := 
      Module[{nullp, p, k, nullprob, prob}, 

        nullp = Exp[nullη]/\((1 + Exp[nullη])\); 

        p = Exp[η]/\((1 + Exp[η])\); 

        nullprob = 
          Table[PDF[BinomialDistribution[m[\([k]\)], nullp[\([k]\)]], 
              y[\([k]\)]], {k, 1, n}]; 

        prob = Table[
            PDF[BinomialDistribution[m[\([k]\)], p[\([k]\)]], y[\([k]\)]], {
              k, 1, n}]; 
{nullprob, prob}
]\)

POISSON

\(poisson[y_, nullη_, η_, n_] := 
      Module[{nullλ, λ, k, nullprob, prob}, 		

        nullλ = Exp[nullη]; 
λ = Exp[η]; 

        nullprob = 
          Table[PDF[PoissonDistribution[nullλ[\([k]\)]], 
              y[\([k]\)]], {k, 1, n}]; 

        prob = Table[
            PDF[PoissonDistribution[λ[\([k]\)]], y[\([k]\)]], {k, 1, 
              n}]; 
{nullprob, prob}
		]\)

GAMMA

\(gamma[y_, nullη_, η_, n_] := 
      Module[{nullm, m, k, nullprob, prob}, 
\ \ \ \ nullm = \ 1/nullη; 
        
		m\  = \ 1/η; 

        nullprob = 
          Table[PDF[GammaDistribution[a, \ nullm[\([k]\)]/a], y[\([k]\)]], {
              k, 1, n}]; 

        prob = Table[
            PDF[GammaDistribution[a, \ m[\([k]\)]/a], y[\([k]\)]], {k, 1, 
              n}]; 
{nullprob, prob}
		]\)

Maximum likelihood estimates

This function returns the exact maximum likelihood estimates of the parameters and their standard deviations

\(MLEst[parlist_, prob_, k2_, family_, opt___] := 
      Module[{nulllogl, logl, nullmaxlogl, maxlogl, nullmle, mle, 
          information, variance, standarderror, i, j}, 
	\ 
        nulllogl = Apply[Plus, Log[prob[\([1]\)]]]; 
	\ 
        logl = Apply[Plus, Log[prob[\([2]\)]]]; 		
\ \ 	
        nullmaxlogl = 
          If[family == "gamma", 
            FindMinimum[\(-nulllogl\), {a, 1}, {b0, 0.0001}], 
            FindMinimum[\(-nulllogl\) /. σ -> 1, {b0, 0.0001}]]; 		

        	\ maxlogl = 
          If[family == "gamma", FindMinimum[\(-logl\), {a, 1}, opt], \ 
            FindMinimum[\(-logl\) /. σ -> 1, opt]]; 		
\ \ \ \ \ 
        nullmle = nullmaxlogl[\([2]\)]; 
		mle = maxlogl[\([2]\)]; 

        information = 
          Table[Table[\ \ 
				
                \(-D[logl, parlist[\([i]\)], parlist[\([j]\)]]\), {j, 1, 
                  k2}], {i, 1, k2}] /. mle; 		

        variance = Inverse[information]; 

        standarderror = Table[Sqrt[variance[\([i, i]\)]], {i, 1, k2}]; 
        
{nullmle, mle, \ standarderror}	
]\)

Deviance

This function returns the deviance residuals and the deviances of  both the specified model and the null model

\(\(	Deviance[y_, nullη_, \ η_, nullparameter_, parameter_, 
        family_] := 
      Module[{nullηhat, ηhat, nullµhat, µhat\ , 
          nulldevres, devres, nulldeviance, deviance}, 

        nullηhat = nullη /. nullparameter; 
		
        ηhat = η /. parameter; 
		
        nullµhat = \ 
          Which[family == "normal", nullηhat, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "binomial", 
            Exp[nullηhat]/\((1 + Exp[nullηhat])\), 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "poisson", Exp[nullηhat], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "gamma", 1/nullηhat]; 	
	
        µhat = \ 
          Which[family == "normal", ηhat, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "binomial", 
            Exp[ηhat]/\((1 + Exp[ηhat])\), 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "poisson", Exp[ηhat], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "gamma", 1/ηhat]; 
		
        nulldevres = 
          Which[family == "normal", 
            Sign[y - nullµhat]\ \((y - nullµhat)\), 

            				\ \ \ \ \ family == "binomial", \ 	
            If[m == Table[1, {Length[m]}], 
              Sign[y - nullµhat]\ 
                \((2 \((
                        \(-y\)\ 
                            Log[\ nullµhat/\((1 - nullµhat)\)]
                               - \ Log[1 - nullµhat\ ])\))\)^0.5, 
              Sign[y - nullµhat\ m]\ 
                \((2 \((
                        y\ Log[\ y/\((nullµhat\ m)\)] + 
                          \((m - y)\)\ 
                            Log[\((m - y)\)/\((m - nullµhat\ m)\)])\))
                    \)^0.5], 
				\ \ \ \ \ family == "poisson", \ 
            Sign[y - nullµhat]\ 
              \((2 \((y\ Log[y/nullµhat] - \((y - nullµhat\ )\))
                      \))\)^0.5, 
				\ \ \ \ \ family == "gamma", 
            \ \ Sign[y - nullµhat]\ 
              \((2 \((\(-Log[\ y/nullµhat]\) + 
                        \((y - nullµhat)\)/nullµhat)\))\)^0.5

            				]; 

        devres = 
          Which[family == "normal", \ \((y - µhat\ )\), 

            				\ \ \ \ \ family == "binomial", \ 		
            If[m == Table[1, {Length[m]}], 
              Sign[y - µhat\ ]\ 
                \((2 \((
                        \(-y\)\ Log[\ µhat/\((1 - µhat)\)\ ] - 
                          Log[1 - µhat\ ])\))\)^0.5, 
				
              Sign[y - µhat\ m] 
                \((2 \((
                        y\ Log[\ y/\((µhat\ m)\)] + 
                          \((m - y)\)\ 
                            Log[\((m - y)\)/\((m - µhat\ m)\)])\))
                    \)^0.5], 
			\ \ family == "poisson", \ 
            Sign[y - µhat]\ 
              \((2 \((y\ Log[y/µhat\ ] - \((y - µhat\ )\))\))
                  \)^0.5, 
				\ \ \ \ \ family == "gamma", \ \ 
            Sign[y - µhat\ ]\ 
              \((2 \((\(-Log[\ y/µhat\ ]\) + 
                        \((y - µhat\ )\)/µhat)\))\)^0.5
		]; 
        		
	nulldeviance = \ Apply[Plus, nulldevres^2]; 
	\ \ 
        deviance = \ Apply[Plus, devres^2]; 
        
{nulldevres, nulldeviance, devres, \ deviance, µhat}
	]
        \)\)

Plot

Diagnositic Plots 
1: plot: Residuals vs Fitted Values
2: plot: Sqrt of abs(Residuals) vs Predictions
3: plot: Response vs Fitted Values
4: plot: Normal QQplot of Std. Residual

Define  Quartile Line for the normal probability plot

\(quartileLine[x_] := 
	
      Module[{slope, \ intercept, \ quartilesX, \ quartilesY, \ extremesX}, 

        		quartilesY = {Quantile[x, 0.25], Quantile[x, 0.75]}; 
		
        quartilesX
           = {Quantile[NormalDistribution[0, 1], 0.25], 
            Quantile[NormalDistribution[0, 1], 0.75]}; 
		
        slope = \((quartilesY[\([2]\)] - quartilesY[\([1]\)])\)/
            \((quartilesX[\([2]\)] - quartilesX[\([1]\)])\); 
		
        intercept = quartilesY[\([1]\)] - slope*quartilesX[\([1]\)]; 

        \ \ \ \ \ 
        extremesX
           = {Quantile[NormalDistribution[0, 1], 0.5/Length[x]], \ \ \ \ \ 
            Quantile[NormalDistribution[0, 1], 1 - 0.5/Length[x]]}; 

        \ \ \ \ \ 
        Line[{{extremesX[\([1]\)], \ intercept + slope*extremesX[\([1]\)]}, 

            		{extremesX[\([2]\)], \ intercept + slope*extremesX[\([2]\)]}}]
          
		]\)

Function to do normal probability plots of a vector of data

\(Nplot[data_] := 
      Module[{p}, p = \((Range[Length[data]] - 0.5)\)/Length[data]; 
        ListPlot[
          Transpose[{Quantile[NormalDistribution[0, 1
				], p], 
              Sort[data]}], PlotRange -> All, 
			
          PlotLabel -> 
            StyleForm["Normal QQplot of Std. Residuals", \ 
              FontSize -> 13, \ FontWeight -> "Bold"], 	
          FrameLabel
             -> {"Quantiles of Standard Normal", 
              "Deviance Residuals"}, \ PlotStyle -> PointSize[0.02], \ 
          Frame -> True, \ Axes -> False, 
			
          Epilog -> quartileLine[data]]]\)

\(DiagPlot[devresiduals_, predict_, fitted_, \ response_]\  := 
      Module[{}, 
\ \ Print["Diagnostic Plots :"]; 		
\ \ 
        ListPlot[Transpose[{fitted, devresiduals}], 	
          FrameLabel -> {"Fitted", "Residuals"}, 
          PlotLabel -> 
            StyleForm["Residuals vs Fitted Values", \ FontSize -> 13, \ 
              FontWeight -> "Bold"], \ PlotStyle -> PointSize[0.02], 
          Frame -> True, \ Axes -> False, 
          Epilog -> {Line[{{Min[fitted] - 2, 0}, {Max[fitted] + 2, 0}}]}

          \ ]; 		
\ 
        ListPlot[Transpose[{predict, Sqrt[Abs[devresiduals]]}], 	
          FrameLabel
             -> {"Predicted", \*""\!\(\@\(|residuals | \)\)""}, 

          			PlotLabel -> 
            StyleForm["Sqrt of |residuals| vs Predictions", \ 
              FontSize -> 13, \ FontWeight -> "Bold"], Frame -> True, \ 
          Axes -> False, PlotStyle -> PointSize[0.02]]; 
	
        If[family != "binomial", m = Table[1, {Length[response]}]]; 
	
        ListPlot[Transpose[{fitted, response/m}], 
          PlotLabel -> 
            StyleForm["Response vs Fitted Values", \ FontSize -> 13, \ 
              FontWeight -> "Bold"], 	
          FrameLabel -> {"Fitted", "Response"}, \ 
          PlotStyle -> PointSize[0.02], \ Frame -> True, \ Axes -> False, 

          			PlotRange
             -> {{Min[fitted, \ response/m] - 0.2, 
                Max[fitted, response/m] + 0.2}, {
                Min[fitted, \ response/m] - 0.2, 
                Max[fitted, response/m] + 0.2}}, 
			
          Epilog -> {
              Line[{{Min[fitted, response/m] - 1, 
                    Min[fitted, response/m] - 1}, {
                    Max[fitted, response/m] + 1, 
                    Max[fitted, response/m] + 1}}]}
\ \ ]; 

        Nplot[devresiduals]
		]\)

ExactGLM

Main Program :
ExactGLM[data,family,fun,var] fits a Generalized linear model of family family to a list of data as a linear combination of the functions funs of 
  variables vars. Canonical links are used as the link functionsfasd

\(ExactGLM[data_, family_, funs_, vars_, plot_] := 
      Module[{k1, k2, n, slots, rules, PureFunction, XMatrixDesignMatrix, y, 
          par = Sequence[], parlist, prob, pos = 0, i\ , newfuns, nullmle, 
          mle, \ standarderror, predict, nulldevres, nulldeviance, devres, \ 
          deviance, µhat, gof, modelchisq, rsq, adjrsq, resulttable, 
          report}, \ 
		k1 = Length[data[\([1]\)]]; 
\ \ \ \ \ 
        y = \(Transpose[data]\)[\([k1]\)]; 
		n = Length[y]; 
	\ \ 
        If[family == "binomial", {m = \(Transpose[data]\)[\([k1 - 1]\)], 
            \ pos = 1}]; 
\ 	\ 
        slots = Table[Slot[n], {n, 1, Length[vars]}]; 
	\ \ 
        rules = Thread[slots -> vars]; 
		
        PureFunction = Map[Release, Function[vars, funs]]; 
\ \ \ \ \ 
        PureFunction = PureFunction /. rules; 
		
        XMatrix = Transpose[Drop[Transpose[data], \(-1\) - pos]]; 
		
        DesignMatrix = Apply[PureFunction, XMatrix, 1]; 

        	\ \ {DesignMatrix2, newfuns} = MakeDummy[DesignMatrix, funs]; 

        				k2 = Length[DesignMatrix2[\([1]\)]]; 
		
        parlist = 
          Table[Symbol[StringJoin["b", ToString[i]]], {i, 0, k2 - 1}]; 

        		nullη = Table[b0, {n}]; 
	\ \ 
        η = DesignMatrix2.parlist; 
		
        prob\  = 
          Which[family == "normal", normal[y, nullη, η, n], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "binomial", binomial[y, nullη, η, m, n], 
            
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "poisson", poisson[y, nullη, η, n], 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            family == "gamma", gamma[y, nullη, η, n]]; 
		
        Do[par = Sequence[par, {parlist[\([i]\)], 0.0001}], {i, 1, k2}]; 

        \ \ \ \ {nullmle, mle, \ standarderror}	 = 
          MLEst[parlist, prob, k2, family, \ par]; 
		
        predict = η /. mle; 

        	{nulldevres, nulldeviance, devres, \ deviance, µhat} = 
          Deviance[y, \ nullη, η, nullmle, mle, family]\ ; 
	
        gof = 1 - CDF[ChiSquareDistribution[n - k2], deviance]; 
	
        modelchisq = 
          1 - CDF[ChiSquareDistribution[k2 - 1], nulldeviance - deviance	]; 
        
		rsq = 1 - deviance/nulldeviance; 
		
        adjrsq = rsq - k2/nulldeviance; 
				
        resulttable = 
          Transpose[{parlist /. mle, 
              standarderror /. σ -> \((deviance/\((n - k2)\))\)^0.5}]; 
        
\ \ \ \ \ 
        report = {DevianceResiduals -> devres, Predict -> predict, \ 
            Fitted -> µhat, 
				
            ParameterTable -> 
              \((TableForm[resulttable, 
                  TableHeadings
                     -> {newfuns, {"Estimate", "Std. Error"}}])\), 

            				ResidualDeviance -> 
              StringJoin[ToString[deviance], " on ", ToString[n - k2], 
                " degrees of freedom"], 
            NullDeviance -> 
              StringJoin[ToString[nulldeviance], " on ", ToString[n - 1], 
                " degrees of freedom"], 
				
            ModelChiSq
               -> {ChiSq -> nulldeviance - deviance, df -> k2 - 1, 
                pvalue -> modelchisq}, 
				RSquared -> rsq, \ 
            AdjustedRSquared -> adjrsq, 
				
            LikelihoodRatioTest -> {df -> n - k2, pvalue -> gof}
		}\ ; 

        		If[plot, DiagPlot[devres, predict, µhat, \ y]]; 
		
        report\ \ \ 
		]\)

Epilogue

End[]