Splus Log related to Fisher-Scoring and Newton-Raphson 
Iteration in a Binomial Regression Setting 
with General Link Function                     4/4/03
======================================================

> binRdat <- rbinom(5, rpois(5,10) -> nvec, 
      0.3 + 0.4*(xvec <- runif(5)))

### These are the data we work with:
> cbind(binRdat, nvec, xvec)
     binRdat nvec      xvec 
[1,]       2   10 0.1289742
[2,]       6   10 0.6965222
[3,]       7   13 0.3797844
[4,]      12   16 0.7889374
[5,]       6   13 0.1761357 

## Note that the actual link function used in
##   generating these data is : identity

> exmpglm1 <- glm(cbind(binRdat, nvec-binRdat) ~ xvec,
    family=binomial) ### default logistic link
Coefficients:
 (Intercept)     xvec 
  -0.9810702 2.487566

Degrees of Freedom: 5 Total; 3 Residual
Residual Deviance: 1.883808 

> exmpglm2 <- glm(cbind(binRdat, nvec-binRdat) ~ xvec,
    family=binomial(link=probit)) ### probit link
Coefficients:
 (Intercept)     xvec 
  -0.6113337 1.544147

Degrees of Freedom: 5 Total; 3 Residual
Residual Deviance: 1.876371 

> exmpglm3B <- glm(binRdat/nvec ~ xvec,
    family=quasi(link=identity, variance= mu),
    weights=nvec) ### linear regression with var=mean
Coefficients:
 (Intercept)     xvec 
   0.2539599 0.613931

Degrees of Freedom: 5 Total; 3 Residual
Residual Deviance: 1.120925 

> exmpglm3 <- glm(binRdat/nvec ~ xvec,
    family=quasi(link=identity, variance= "mu(1-mu)"),
    weights=nvec) ### identity link: binomial form
### specified by the choice of variance as fcn of mu
Coefficients:
 (Intercept)      xvec 
   0.2608657 0.5955406

Degrees of Freedom: 5 Total; 3 Residual
Residual Deviance: 1.846349 


### Now we code an Splus function to take data from this 
###   structure with general inverse-link-function $G$ and 
###   initial guess for beta parameters and do a 
###   Newton-Raphson and Fisher Scoring iteration.

> MLhistBin <- function(binit, size, resp, prdctr, 
     G = plogis, Gpr = dlogis, Gdpr = function(x) 
     dlogis(x)*(1-2*plogis(x)))
{
eta <- binit[1]+prdctr*binit[2]
mu <- G(eta)
aux <- (resp-size*mu)*Gpr(eta)/(mu*(1-mu))
scor <- c(sum(aux), sum(prdctr*aux))
aux2 <- size*(Gpr(eta)^2)/(mu*(1-mu))
FInfo <- matrix(c(sum(aux2),rep(sum(aux2*prdctr),2),
     sum(aux2*prdctr^2)), ncol=2)
aux3 <- ((resp - size * mu) * ((1-2*mu)*Gpr(eta)^2 - 
     (mu-mu^2)*Gdpr(eta))/(mu*(1-mu))^2
Info <- FInfo + matrix(c(sum(aux3), rep(sum(aux3*
     prdctr),2), sum(aux3*prdctr^2)), ncol=2)
list(BnewNR = binit + solve(Info, scor), BnewFS = binit +
     solve(FInfo, scor), Bold = binit)
}

> nxt <- MLhistBin(rep(0,2), nvec, binRdat, xvec)
$BnewNR:
[1] -0.9394949  2.3571914
$BnewFS:
[1] -0.9394949  2.3571914
$Bold:
[1] 0 0

> nxtA <- nxt$BnewNR
> for(i in 1:5) { 
   nxtA <- MLhistBin(nxtA, nvec, binRdat, xvec)$BnewNR
   cat( round(nxtA,4), "\n") }
-0.9807 2.4863 
-0.9811 2.4876 
-0.9811 2.4876 
-0.9811 2.4876 
-0.9811 2.4876 
> nxtB <- nxt$BnewFS
> for(i in 1:5) { 
   nxtB <- MLhistBin(nxtB, nvec, binRdat, xvec)$BnewFS
   cat( round(nxtB,4), "\n") }
-0.9807 2.4863 
-0.9811 2.4876 
-0.9811 2.4876 
-0.9811 2.4876 
-0.9811 2.4876 

> nxt <- MLhistBin(rep(0,2), nvec, binRdat, xvec, G=pnorm,
    Gpr = dnorm, Gdpr = function(x) -x*dnorm(x))
$BnewNR:
[1] -0.5887411  1.4771506
$BnewFS:
[1] -0.5887411  1.4771506
$Bold:
[1] 0 0

> nxtA <- nxt$BnewNR
  for(i in 1:5) { 
   nxtA <- MLhistBin(nxtA, nvec, binRdat, xvec, G=pnorm,
     Gpr = dnorm, Gdpr = function(x) -x*dnorm(x))$BnewNR
   cat( round(nxtA,4), "\n") }
-0.6113 1.5439 
-0.6113 1.5441 
-0.6113 1.5441 
-0.6113 1.5441 
-0.6113 1.5441
> nxtB <- nxt$BnewFS
  for(i in 1:5) { 
   nxtB <- MLhistBin(nxtB, nvec, binRdat, xvec, G=pnorm,
     Gpr = dnorm, Gdpr = function(x) -x*dnorm(x))$BnewFS
   cat( round(nxtB,4), "\n") }
-0.6111 1.5435 
-0.6113 1.5441 
-0.6113 1.5441 
-0.6113 1.5441 
-0.6113 1.5441

> nxt <- MLhistBin(rep(0.5,2), nvec, binRdat, xvec, G= 
    function(x) x , Gpr = function(x) 1, Gdpr = function(x) 0)
$BnewNR:
[1] 0.2742102 0.6978773
$BnewFS:
[1] 0.2645638 0.5891888
$Bold:
[1] 0.5 0.5

> nxtA <- nxt$BnewNR
  for(i in 1:5) { 
   nxtA <- MLhistBin(nxtA, nvec, binRdat, xvec, G= function(x) 
    x , Gpr = function(x) 1, Gdpr = function(x) 0)$BnewNR
   cat( round(nxtA,4), "\n") }
0.2501 0.6411 
0.2603 0.5979 
0.2609 0.5955 
0.2609 0.5955 
0.2609 0.5955 
> nxtB <- nxt$BnewFS
  for(i in 1:5) { 
   nxtB <- MLhistBin(nxtB, nvec, binRdat, xvec, G= function(x) 
    x , Gpr = function(x) 1, Gdpr = function(x) 0)$BnewFS   
   cat( round(nxtB,4), "\n") }
0.261 0.5953 
0.2609 0.5955 
0.2609 0.5955 
0.2609 0.5955 
0.2609 0.5955

### NOTE that the Splus solutions in every case agree with the 
NR and Fisher-score soluttions after just 2 iterations !!

### NOTE also that in the canonical-link case, the logistic,
Fisher scoring and NR agree perfectly in all iterations, 
but the NR and Fisher-score solution differ on the very first
iteration in the non-canonical-link cases.