HW9comp
Gillian Cheng
2023-12-05
FYI
Ex 20-4:
- Please use the Newton-Raphson method to find the maximum likelihood estimate (MLE) of the regression coefficients of poisson regression.
rate <- read.csv("rate.csv")
rate$Age.f <- factor(rate$Age)
# define X and Y
Y <- rate$Death # 1 by 24
dummy <- diag(1, length(levels(rate$Age.f))-1)
dummy <- rbind(rep(0,length(levels(rate$Age.f))-1),dummy)
dummy <- rbind(dummy, dummy)
X <- cbind(rep(1,length(Y)),dummy, ifelse(rate$sex == "m",1,0)) #24 by 13
ftn3 <- function(betacoef){
mu <- exp(X%*%betacoef+log(rate$PY/100000))
gradient <- t(X)%*%(Y-mu)
Hessian <- -t(X)%*%diag(c(mu),length(Y))%*%X
logL <- sum(-mu + Y*log(mu)-log(factorial(Y)))
return(list(gradient, Hessian, logL))
}(beta3 <- newtonraphson(ftn3,c(0,0,0,0,0,0,0,0,0,0,0,0,0)))## Algorithm converged## [,1]
## [1,] 0.94453534
## [2,] -0.11158760
## [3,] -0.19307396
## [4,] -0.39994535
## [5,] -0.57528753
## [6,] -0.40114660
## [7,] -0.31142104
## [8,] 0.04452828
## [9,] 0.07301330
## [10,] 0.17111769
## [11,] 0.15256371
## [12,] -0.27115727
## [13,] 0.56064525cc <- glm(Death ~ Age.f + sex, offset = log(PY/100000), data = rate, family = "poisson")
cc$coefficients## (Intercept) Age.f2 Age.f3 Age.f4 Age.f5 Age.f6
## 0.94453534 -0.11158760 -0.19307396 -0.39994535 -0.57528753 -0.40114660
## Age.f7 Age.f8 Age.f9 Age.f10 Age.f11 Age.f12
## -0.31142104 0.04452828 0.07301330 0.17111769 0.15256371 -0.27115727
## sexm
## 0.560645252.Please find the variance-covariance matrix for the betas .
solve(-ftn3(beta3)[[2]])## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.014509599 -1.084465e-02 -1.086946e-02 -1.089294e-02 -1.090877e-02
## [2,] -0.010844653 2.369025e-02 1.086957e-02 1.086941e-02 1.086930e-02
## [3,] -0.010869462 1.086957e-02 2.726301e-02 1.086956e-02 1.086956e-02
## [4,] -0.010892943 1.086941e-02 1.086956e-02 3.718550e-02 1.086982e-02
## [5,] -0.010908767 1.086930e-02 1.086956e-02 1.086982e-02 4.933153e-02
## [6,] -0.010942335 1.086907e-02 1.086956e-02 1.087003e-02 1.087035e-02
## [7,] -0.011051674 1.086832e-02 1.086956e-02 1.087073e-02 1.087153e-02
## [8,] -0.011099002 1.086799e-02 1.086956e-02 1.087104e-02 1.087204e-02
## [9,] -0.011212934 1.086722e-02 1.086956e-02 1.087177e-02 1.087326e-02
## [10,] -0.011420939 1.086579e-02 1.086955e-02 1.087311e-02 1.087550e-02
## [11,] -0.011707698 1.086383e-02 1.086954e-02 1.087495e-02 1.087859e-02
## [12,] -0.012212644 1.086037e-02 1.086953e-02 1.087819e-02 1.088403e-02
## [13,] -0.005718767 -3.913833e-05 -1.619563e-07 3.672866e-05 6.158833e-05
## [,6] [,7] [,8] [,9] [,10]
## [1,] -0.0109423347 -0.0110516743 -0.0110990021 -0.0112129343 -0.0114209388
## [2,] 0.0108690672 0.0108683189 0.0108679950 0.0108672153 0.0108657917
## [3,] 0.0108695632 0.0108695601 0.0108695587 0.0108695555 0.0108695496
## [4,] 0.0108700326 0.0108707348 0.0108710388 0.0108717705 0.0108731064
## [5,] 0.0108703489 0.0108715264 0.0108720361 0.0108732631 0.0108755032
## [6,] 0.0453537786 0.0108732058 0.0108741520 0.0108764297 0.0108805880
## [7,] 0.0108732058 0.0431367406 0.0108810438 0.0108867438 0.0108971502
## [8,] 0.0108741520 0.0108810438 0.0358840270 0.0108912083 0.0109043191
## [9,] 0.0108764297 0.0108867438 0.0108912083 0.0431600202 0.0109215770
## [10,] 0.0108805880 0.0108971502 0.0109043191 0.0109215770 0.0509530845
## [11,] 0.0108863207 0.0109114966 0.0109223940 0.0109486274 0.0109965213
## [12,] 0.0108964153 0.0109367588 0.0109542215 0.0109962596 0.0110730079
## [13,] 0.0001143263 0.0002861071 0.0003604625 0.0005394587 0.0008662494
## [,11] [,12] [,13]
## [1,] -0.01170770 -0.012212644 -5.718767e-03
## [2,] 0.01086383 0.010860373 -3.913833e-05
## [3,] 0.01086954 0.010869527 -1.619563e-07
## [4,] 0.01087495 0.010878191 3.672866e-05
## [5,] 0.01087859 0.010884030 6.158833e-05
## [6,] 0.01088632 0.010896415 1.143263e-04
## [7,] 0.01091150 0.010936759 2.861071e-04
## [8,] 0.01092239 0.010954221 3.604625e-04
## [9,] 0.01094863 0.010996260 5.394587e-04
## [10,] 0.01099652 0.011073008 8.662494e-04
## [11,] 0.07772922 0.011178815 1.316770e-03
## [12,] 0.01117881 0.154222269 2.110077e-03
## [13,] 0.00131677 0.002110077 8.984612e-03vcov(cc)## (Intercept) Age.f2 Age.f3 Age.f4
## (Intercept) 0.014509599 -1.084465e-02 -1.086946e-02 -1.089294e-02
## Age.f2 -0.010844653 2.369025e-02 1.086957e-02 1.086941e-02
## Age.f3 -0.010869462 1.086957e-02 2.726301e-02 1.086956e-02
## Age.f4 -0.010892943 1.086941e-02 1.086956e-02 3.718550e-02
## Age.f5 -0.010908767 1.086930e-02 1.086956e-02 1.086982e-02
## Age.f6 -0.010942335 1.086907e-02 1.086956e-02 1.087003e-02
## Age.f7 -0.011051674 1.086832e-02 1.086956e-02 1.087073e-02
## Age.f8 -0.011099002 1.086799e-02 1.086956e-02 1.087104e-02
## Age.f9 -0.011212934 1.086722e-02 1.086956e-02 1.087177e-02
## Age.f10 -0.011420939 1.086579e-02 1.086955e-02 1.087311e-02
## Age.f11 -0.011707698 1.086383e-02 1.086954e-02 1.087495e-02
## Age.f12 -0.012212644 1.086037e-02 1.086953e-02 1.087819e-02
## sexm -0.005718766 -3.913832e-05 -1.619563e-07 3.672866e-05
## Age.f5 Age.f6 Age.f7 Age.f8
## (Intercept) -1.090877e-02 -0.0109423347 -0.0110516743 -0.0110990020
## Age.f2 1.086930e-02 0.0108690672 0.0108683189 0.0108679950
## Age.f3 1.086956e-02 0.0108695631 0.0108695601 0.0108695587
## Age.f4 1.086982e-02 0.0108700326 0.0108707348 0.0108710388
## Age.f5 4.933152e-02 0.0108703489 0.0108715264 0.0108720361
## Age.f6 1.087035e-02 0.0453537786 0.0108732058 0.0108741520
## Age.f7 1.087153e-02 0.0108732058 0.0431367391 0.0108810438
## Age.f8 1.087204e-02 0.0108741520 0.0108810438 0.0358840269
## Age.f9 1.087326e-02 0.0108764296 0.0108867438 0.0108912083
## Age.f10 1.087550e-02 0.0108805880 0.0108971502 0.0109043191
## Age.f11 1.087859e-02 0.0108863207 0.0109114966 0.0109223940
## Age.f12 1.088403e-02 0.0108964153 0.0109367588 0.0109542215
## sexm 6.158832e-05 0.0001143263 0.0002861071 0.0003604625
## Age.f9 Age.f10 Age.f11 Age.f12 sexm
## (Intercept) -0.0112129343 -0.0114209388 -0.01170770 -0.012212644 -5.718766e-03
## Age.f2 0.0108672152 0.0108657917 0.01086383 0.010860373 -3.913832e-05
## Age.f3 0.0108695555 0.0108695496 0.01086954 0.010869527 -1.619563e-07
## Age.f4 0.0108717705 0.0108731064 0.01087495 0.010878191 3.672866e-05
## Age.f5 0.0108732631 0.0108755032 0.01087859 0.010884030 6.158832e-05
## Age.f6 0.0108764296 0.0108805880 0.01088632 0.010896415 1.143263e-04
## Age.f7 0.0108867438 0.0108971502 0.01091150 0.010936759 2.861071e-04
## Age.f8 0.0108912083 0.0109043191 0.01092239 0.010954221 3.604625e-04
## Age.f9 0.0431600202 0.0109215770 0.01094863 0.010996260 5.394586e-04
## Age.f10 0.0109215770 0.0509530837 0.01099652 0.011073008 8.662493e-04
## Age.f11 0.0109486274 0.0109965213 0.07772909 0.011178815 1.316770e-03
## Age.f12 0.0109962595 0.0110730078 0.01117881 0.154222064 2.110077e-03
## sexm 0.0005394586 0.0008662493 0.00131677 0.002110077 8.984611e-03- Please find the log likelihood at the betas .
ftn3(beta3)[[3]]## [1] -59.38966logLik(cc)## 'log Lik.' -59.38966 (df=13)