Step-by-step formula based _R_ codes for GLM and GLMM

I know how to fit generalized linear models (GLMs) and generalized linear mixed models (GLMMs) with glm and glmer from lme4 package in R. Being a student of statistics, I'm interested in learning how to fit GLM and GLMM following step-by-step formula bases R codes. I'd highly appreciate if you point开发者_Go百科 out any resource and/or reference in this regard. Thanks in advance.

EDiT

I'd like to do GLM and GLMM step by step using formula as we do LM using matrix approach. Is there any R book or tutorial available that use this approach? Thanks

"An R Companion to Applied Regression" by Fox and Weisberg, has an excellent guide in chapter 8, with logistic regression as an example. The book also teaches a bit about how to create model functions in general with S3 and S4 objects. In particular, it has good answers to a recent question I'd asked about modeling -- What are the key components and functions for standard model objects in R?.

This may help
** Poisson regression: GLM**
*Suggested reading: An Introduction to generalized linear model, by Annette J. Dobson, 2nd Edition, Chapter 4, section 4.3 and 4.4 *

library(MASS)
poisreg = function(n, b1, y, x1, tolerence) {  # n is the number of iteration   
  x0 = rep(1, length(x1))   
  x = cbind(x0, x1)  
  y = as.matrix(y)  
  w = matrix(0, nrow = (length(y)), ncol = (length(y)))  
  b0 = b1  
  result = b0
  for (i in 1:n) {  
    mu = exp(x %*% b0)     
    diag(w) = mu  
    eta = x %*% b0  
    z = eta + (y - mu) * (1/mu)   # dot product of (y - mu) & (1/mu)   
    xtwx = t(x) %*% w %*% x  
    xtwz = t(x) %*% w %*% z  
    b1 = solve(xtwx, xtwz)  
    if(sqrt(sum(b0 - b1)^2) > tolerence) (b0 <- b1)  
    result<- cbind(result,b1) # to get all the iterated values  
  }  
  result  
}
x1 <- c(-1,-1,0,0,0,0,1,1,1) # x1 is the explanatory variable 
y<- c(2,3,6,7,8,9,10,12,15)  # y is the dependent variable
b1 = c(1,2) # initial value  
poisreg (10, b1, y, x1, .001)   # Nicely converge after 10 iterations  
glm(y~x1, family=poisson(link="log"))   # check your result with the R GLM program