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

Question

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 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

Do you mean that you want to learn how to write code to fit GLM(M)s? — Gavin Simpson, Aug 03 '11 at 16:41
I think the answer is, in part, "McCullagh and Nelder". Read that, tells you all the algorithms. Start with simple linear Gaussian stuff first though. — Spacedman, Aug 03 '11 at 16:46
Crawley's "The R book" has an entire chapter on GLM using `glm`and a 5-page worked example on GLMM usig `lmer` — Andrie, Aug 03 '11 at 16:52
This question contains a list of possibly useful resources: http://stackoverflow.com/q/192369/602276 — Andrie, Aug 04 '11 at 07:35
For GLMs Faraway's book on *Extending the linear model* is good. For GLMMs there's not much out there -- see Zuur et al (Mixed models) and http://glmm.wikidot.com/faq — Ben Bolker, Nov 29 '11 at 22:22

overwhelmed · Answer 1 · 2013-02-01T16:40:16.107

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

+1: Thanks @overwhelmed for your nice answer. You answer could be good start. Thanks again. — MYaseen208, Jan 31 '13 at 11:13
Is there a way to change this code for a binomial implementation instead of poisson? — James H, Jan 11 '16 at 03:29

score 2 · Answer 2 · edited May 23 '17 at 12:01

"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?.

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

2 Answers2

Linked