• API
• FAQ
• Tools
• Archive
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.

Top