SHARE
TWEET

Untitled

a guest Jul 21st, 2019 89 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ## generate random data
  2. y <- rpois(25, 2)
  3. x <- rnorm(25) + (log(replace(y, y == 0, 0.25)))
  4.      
  5. mod1 <- optim(c(1, 1), function(p)
  6.     with(list(mu=p[1] + p[2] * x), -sum(y * mu - exp(mu))))
  7.  
  8. mod2 <- glm(y ~ x, family='poisson', data=data.frame(y=y, x=x))
  9.  
  10. ## check that the parameters coincide
  11. mod1$par
  12. coefficients(mod2)
  13.      
  14. ## eta = a + bx
  15. eta <- cbind(1, x) %*% mod1$par
  16. tail(eta)
  17. tail(predict(mod2))
  18.  
  19. ## yhat = exp(eta)
  20. yhat <- exp(eta)
  21. tail(yhat)
  22. tail(fitted.values(mod2))
  23.  
  24. ## log-likelihood
  25. sum(dpois(y, yhat, log=TRUE))
  26. logLik(mod2)[1]
  27.      
  28. ## deviance: D = 2(logLik_fullmodel - logLik_model)
  29. 2 * (sum(dpois(y, y, log=TRUE)) - sum(dpois(y, yhat, log=TRUE)))
  30. mod2$deviance
  31.      
  32. tail(2 * (yhat - y - y * (eta - log(y))))
  33. tail(resid(mod2)^2)
  34. sum(resid(mod2)^2)
  35.      
  36. barplot(rbind(dpois(0:4, y[1], log=TRUE),
  37.               dpois(0:4, yhat[1], log=TRUE)),
  38.         beside=TRUE, names.arg=paste0('y=', 0:4),
  39.         legend.text=c('saturated', 'model'), args.legend=list(x='bottom'))
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top