function slide=glmdemo
% This is a slideshow file for use with playshow.m and makeshow.m
% To see it run, type 'playshow glmdemo', 

% Copyright 1993-2000 The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 2000/05/26 18:52:55 $
if nargout<1,
  playshow glmdemo
else
  %========== Slide 1 ==========

  slide(1).code={
   'subplot(''Position'',[0.10 0.47 0.65 0.45]);',
   'x=1:20;',
   'y = normcdf(2.66 - 0.31 * x);',
   'z = [0.007 0.008 0.067 0.045 0.19 0.088 0.1 0.28 0.65 0.82 0.69 0.95 0.91 0.96 0.97 0.98 1 1 1 0.99];',
   'plot(x,y,''b-'',   ''LineWidth'',2);',
   'hold on;',
   'plot(x,1-z,''r.'', ''MarkerSize'',12);',
   'hold off',
   'text(8,.8,''Prob[success] = \Phi(2.66 - 0.31x)'')',
   'title(''Probit model'');',
   '' };
  slide(1).text={
   'Press the Start button for a demonstration of generalized linear models.',
   '',
   'This demonstration uses the Statistics Toolbox function GLMFIT, which can fit models like the probit model shown here.'};

  %========== Slide 2 ==========

  slide(2).code={
   'subplot(''Position'', [.1 .75 .3 .2]); plot(-1.2:.05:2.5, normcdf(-1.2:.05:2.5),''r-'')',
   'set(gca,''xlim'',[-1.2 2.5])',
   'title(''Logistic function'');',
   '',
   'subplot(''Position'', [.1 .45 .3 .2]); plot(.5:.01:2, 1./(.5:.01:2),''r-'')',
   'title(''Reciprocal function'');',
   '',
   'subplot(''Position'', [.45 .75 .3 .2]); hist(binornd(100,.98,100,1),90:1:100)',
   'title(''Binomial distribution'');',
   '',
   'set(gca,''xlim'',[89.5 100.5])',
   'subplot(''Position'', [.45 .45 .3 .2]);  hist(gamrnd(5,.1,100,1),.1:.2:2)',
   'set(gca,''xlim'',[0 3.2])',
   'title(''Poisson distribution'');' };
  slide(2).text={
   'Ordinary least squares regression fits straight lines, and other functions that are linear in the parameters, to normally distributed data.',
   '',
   'Generalized linear models support other types of functions, known as "link" functions.  Examples of useful link functions include the S-shaped logistic function and the reciprocal function.',
   '',
   'Generalized linear models also support data having various distributions including the binomial, gamma, normal, and Poisson distributions.'};

  %========== Slide 3 ==========

  slide(3).code={
   'subplot(''Position'',[0.10 0.47 0.65 0.45]);',
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   'plot(x, p, ''s'')',
   'xlabel Weight',
   'ylabel Proportion' };
  slide(3).text={
   'This graph shows the proportion of cars of various weights that fail a mileage test.  Our data consist of:',
   '  x = a set of car weights',
   '  n = the number of cars at each weight',
   '  y = the number of cars failing the test at each weight',
   '  p = y ./ n',
   '',
   'The graph is a scatter plot of p as a function of x.',
   ''};

  %========== Slide 4 ==========

  slide(4).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   'b1 = polyfit(x,p,1);',
   'plot(x,p,''s'', x, polyval(b1,x), ''r-'',[2000 4500],[0 0],''k:'',[2000 4500],[1 1],''k:'')',
   'xlabel Weight',
   'ylabel Proportion' };
  slide(4).text={
   'It''s reasonable to assume each y value has a binomial distribution, with a P parameter that depends on x.  But how does P depend on weight?',
   '',
   'We can try fitting a straight line to this plot, but there are two problems:',
   '',
   '1.  The line predicts proportions less than 0 and greater than 1.',
   '2.  The points are clearly not normally distributed, so they violate one of the assumptions usually required for fitting regression models.'};

  %========== Slide 5 ==========

  slide(5).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   'b3 = polyfit(x,p,3);',
   'plot(x,p,''s'', x, polyval(b3,x), ''r-'')',
   'xlabel Weight',
   'ylabel Proportion' };
  slide(5).text={
   'Using a higher-order polynomial may appear to help, but it still has the same problems.  This graph shows that the fitted proportion starts to decrease as Weight goes above 4000; in fact it will become negative for larger Weight values.',
   '',
   'Of course, the assumption of a normal distribution is still incorrect.'};

  %========== Slide 6 ==========

  slide(6).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   'b = glmfit(x,[y n],''binomial'');',
   'yfit = 1./(1+exp(-(b(1) + b(2) * x)));',
   'plot(x,p,''bs'', x,yfit,''r-'')' };
  slide(6).text={
   'Instead, a good approach is to use the GLMFIT function to fit a logistic regression model.  It uses a fitting method that is appropriate for the binomial distribution.  Also, the logistic function is limited to the desirable range [0,1].',
   '',
   '>> [b,d] = glmfit(x,[y n],''binomial'');',
   '>> yfit = 1./(1+exp(-(b(1) + b(2) * x)));',
   '>> plot(x,p,''bs'', x,yfit,''r-'')'};

  %========== Slide 7 ==========

  slide(7).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   '[b1,d1] = glmfit(x,[y n],''binomial'');',
   '[b2,d2] = glmfit([x x.^2], [y n], ''binomial'');',
   'd1 - d2;',
   'pval = 1 - chi2cdf(ans,1);',
   'y1 = 1./(1+exp(-(b1(1) + b1(2) * x)));',
   'y2 = 1./(1+exp(-(b2(1) + b2(2) * x + b2(3) * x.^2)));',
   'plot(x,p,''bs'', x,y1,''r-'', x,y2,''g--'')' };
  slide(7).text={
   'The GLMFIT function provides a number of outputs for examing and testing the model.  For example, we can compare the deviance values for two models to determine if a squared term would improve the fit significantly.  Here, it does not:',
   '',
   '>> [b2,d2] = glmfit([x x.^2], [y n], ''binomial'');',
   '>> diff = d - d2,   pval = 1 - chi2cdf(ans,1)',
   'diff =',
   '    0.7027',
   'pval =',
   '    0.4019'};

  %========== Slide 8 ==========

  slide(8).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   '[b,d,s] = glmfit(x,[y n],''binomial'');',
   'normplot(s.residp)' };
  slide(8).text={
   'To check the model, we can look at a probability plot of the Pearson residuals.  These are normalized to have roughly a standard normal distribution.  (Without this standardization, the residuals would have different variances.)',
   '',
   '>> [b,d,s] = glmfit(x,[y n],''binomial'');',
   '>> normplot(s.residp)',
   '',
   'The normal probability plot does shows a nice agreement with the normal distribution.'};

  %========== Slide 9 ==========

  slide(9).code={
   'x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'';',
   'y = [1 2 0 3 8 8 14 17 19 15 17 21]'';',
   'n = [48 42 31 34 31 21 23 23 21 16 17 21]'';',
   'p = y ./ n;',
   '[b,d,s] = glmfit(x,[y n],''binomial'');',
   'xx = 2500:500:4000;',
   '[yhat,dylo,dyhi] = glmval(b,xx,''logit'',s,.95,100);',
   'errorbar(xx,yhat,dylo,dyhi,'':'');' };
  slide(9).text={
   'Once we are satisfied with the model, we can use it to make predictions.',
   '',
   'Here we estimate the expected number of cars, out of 100 tested, that would pass the mileage test at each of four weights:',
   '',
   '>> xx = 2500:500:4000;',
   '>> [yhat,dylo,dyhi] = glmval(b,xx,''logit'',s,.95,100);',
   '>> errorbar(xx,yhat,dylo,dyhi,'':'');'};
end