Fall 200 STA 216 October 4, 2000Fisher Scoring and Diagnostics1 Fisher ScoringThe Fisher Scoring algorithm can be implemented using weighted least squares regression routines. Given astarting value for β (or π), we construct the vector of weights W and the working response Z, and then findˆβ by regressing Z on X using weights W . We construct new weights and working responses Z, and theniterate until we obtain convergence.For the logit link in binary regression• Z = log(π/(1 − π)) + (Y − π)/(π ∗ (1 − π))• W = π ∗ (1 − π)Here is a no-frills function to perform k iterations of Fisher Scoring under the logit link(equivalent toNewton Raphson because observed and expected information are the same). The input X is assumed tocontain a column of ones. No error checking or comments! (bad prof)Fisher.it _ function(Y,X, pi0, niter=1, print=F) {pi _ pi0for (i in 1:niter) {W _ pi*(1-pi)Z _ log(pi/(1-pi)) + (Y - pi)/(pi*(1-pi))lmobj _ lm(Z ~ X - 1, weights=W)beta _ lmobj$coefeta _ X %*% betapi _ exp(eta)/(1 + exp(eta))if (print) {print(paste("Iteration ", as.character(i), ": Betahat"))print(beta)}}XWX _ t(lmobj$R) %*% lmobj$Rreturn(beta, XWX, pi, W)}Start with initial probabilities pi0 = 0.75 if Y=1, 0.25 if Y=0.After one iteration:> out _ Fisher.it(orings$failure, X, pi0, 1)> out$betaX1 X29.422777 -0.14926471Run for 4 iterations:> out _ Fisher.it(orings$failure, X, pi0, 4, print=T)[1] "Iteration 1 : Betahat"X1 X29.422777 -0.1492647[1] "Iteration 2 : Betahat"X1 X210.76226 -0.1696113[1] "Iteration 3 : Betahat"X1 X210.87462 -0.1713095[1] "Iteration 4 : Betahat"X1 X210.87535 -0.1713205># Asymptotic Covariance Matrix> solve(out$XWX)X1 X2X1 32.5232111 -0.473880521X2 -0.4738805 0.006962009# Standard errors> sqrt(diag(solve(out$XWX)))[1] 5.70291251 0.08343865# correlation between betahat0 and betahat1> -.473880521/( 5.70291251*0.08343865)[1] -0.9958751#Splus:> summary(oring.logit)Coefficients:Value Std. Error t value(Intercept) 10.8753321 5.69793801 1.908643temp -0.1713202 0.08336339 -2.055102Number of Fisher Scoring Iterations: 4Correlation of Coefficients:(Intercept)temp -0.9958713Pretty close! difference in SE(beta0)?22 Residuals and Diagnostics• ordinary residuals(yi− µi)• Person residuals((yi− µi)/pV ar(yi))• Deviance residuals ( sum squared deviance residuals = deviance)sign(˜ηi− ˆηi)p2l(˜ηi, yi) − 2l(ˆηi, yi)• Cook’s Distance(ˆβ(i)−ˆβ)0X0W X(ˆβ(i)−ˆβ)/(p) ≈ (1/p)(eF)2hii/(1 − hii)Studentized forms of the Pearson and Deviance residuals are defined as eFand eDby taking the above anddividing by the leveragehii= diag(H) H = W1/2X(X0W X)−1X0W1/2.In Splus, the residuals function will return each of the un-standardized residuals. Cook’s Distance canbe approximated using the linear model results applied at the last step of the Fisher Scoring algorithm. Thefunction below uses the lm.influence function, which provides estimates of β with single case deletions,and the leverage values.> CooksD _ function(fit){fit.s <- summary.glm(fit)fit.infl <- lm.influence(fit)R <- fit$RXWX <- t(R) %*% RD <- rep(0, length(fit$residuals))Dcoef <- scale(fit.infl$coefficients, center = fit$coefficients, scale= F)disper <- fit.s$dispersion#estimate of stdev^2 for Gaussianp <- sum(fit.infl$hat)for(i in 1:length(D)) {D[i] <- (t(Dcoef[i, ]) %*% XWX %*% Dcoef[i, ])/(p * disper)}D}Construct diagnostics and plots:> h _ lm.influence(oring.logit)$hat> pearson.residual _ residuals(oring.logit, type="pearson")/sqrt(1-h)> dev.residual _ residuals(oring.logit, type="deviance")/sqrt(1-h)> cooksd _ CooksD(oring.logit)> par(mfrow=c(2,2))> index _ 1:24> plot(index, h)> qqnorm(pearson.residual)> plot(index, cooksd)> qqnorm(dev.residual)3••••••••••••••••••••••••indexh5 10 200.05 0.10 0.15 0.20••••••••••••• •••••••••••Quantiles of Standard Normaldev.residual-2 -1 0 1 2-1 0 1 2••••••••••••••••••••••••indexcooksd5 10 200.0 0.10 0.20••••••••••••• •••••••••••Quantiles of Standard Normaldev.residual-2 -1 0 1 2-1 0 1 2Diagnostics for the oring data using the logit link.4•••••••••••••••••••fitted(oring.logit)pearson.residual0.2 0.4 0.6 0.8-1 0 1 2Studentized Pearson Residuals•••••••••••••••••••fitted(oring.logit)dev.residual0.2 0.4 0.6 0.8-1 0 1 2Studentized Deviance ResidualsResidual plots for the oring data using the logit