<<data>>
id	y	visit	trt	bline	age	x1	ltime
101	76	0	1	76	18	0	2.079441542
101	11	1	1	76	18	1	0.693147181
101	14	2	1	76	18	1	0.693147181
101	9	3	1	76	18	1	0.693147181
101	8	4	1	76	18	1	0.693147181
102	38	0	1	38	32	0	2.079441542
102	8	1	1	38	32	1	0.693147181
102	7	2	1	38	32	1	0.693147181
102	9	3	1	38	32	1	0.693147181
102	4	4	1	38	32	1	0.693147181
103	19	0	1	19	20	0	2.079441542
103	0	1	1	19	20	1	0.693147181
103	4	2	1	19	20	1	0.693147181
103	3	3	1	19	20	1	0.693147181
103	0	4	1	19	20	1	0.693147181
104	11	0	0	11	31	0	2.079441542
104	5	1	0	11	31	1	0.693147181
104	3	2	0	11	31	1	0.693147181
104	3	3	0	11	31	1	0.693147181
104	3	4	0	11	31	1	0.693147181
106	11	0	0	11	30	0	2.079441542
106	3	1	0	11	30	1	0.693147181
106	5	2	0	11	30	1	0.693147181
106	3	3	0	11	30	1	0.693147181
106	3	4	0	11	30	1	0.693147181
107	6	0	0	6	25	0	2.079441542
107	2	1	0	6	25	1	0.693147181
107	4	2	0	6	25	1	0.693147181
107	0	3	0	6	25	1	0.693147181
107	5	4	0	6	25	1	0.693147181
108	10	0	1	10	30	0	2.079441542
108	3	1	1	10	30	1	0.693147181
108	6	2	1	10	30	1	0.693147181
108	1	3	1	10	30	1	0.693147181
108	3	4	1	10	30	1	0.693147181
110	19	0	1	19	18	0	2.079441542
110	2	1	1	19	18	1	0.693147181
110	6	2	1	19	18	1	0.693147181
110	7	3	1	19	18	1	0.693147181
110	4	4	1	19	18	1	0.693147181
111	24	0	1	24	24	0	2.079441542
111	4	1	1	24	24	1	0.693147181
111	3	2	1	24	24	1	0.693147181
111	1	3	1	24	24	1	0.693147181
111	3	4	1	24	24	1	0.693147181
112	31	0	1	31	30	0	2.079441542
112	22	1	1	31	30	1	0.693147181
112	17	2	1	31	30	1	0.693147181
112	19	3	1	31	30	1	0.693147181
112	16	4	1	31	30	1	0.693147181
113	14	0	1	14	35	0	2.079441542
113	5	1	1	14	35	1	0.693147181
113	4	2	1	14	35	1	0.693147181
113	7	3	1	14	35	1	0.693147181
113	4	4	1	14	35	1	0.693147181
114	8	0	0	8	36	0	2.079441542
114	4	1	0	8	36	1	0.693147181
114	4	2	0	8	36	1	0.693147181
114	1	3	0	8	36	1	0.693147181
114	4	4	0	8	36	1	0.693147181
116	66	0	0	66	22	0	2.079441542
116	7	1	0	66	22	1	0.693147181
116	18	2	0	66	22	1	0.693147181
116	9	3	0	66	22	1	0.693147181
116	21	4	0	66	22	1	0.693147181
117	11	0	1	11	27	0	2.079441542
117	2	1	1	11	27	1	0.693147181
117	4	2	1	11	27	1	0.693147181
117	0	3	1	11	27	1	0.693147181
117	4	4	1	11	27	1	0.693147181
118	27	0	0	27	29	0	2.079441542
118	5	1	0	27	29	1	0.693147181
118	2	2	0	27	29	1	0.693147181
118	8	3	0	27	29	1	0.693147181
118	7	4	0	27	29	1	0.693147181
121	67	0	1	67	20	0	2.079441542
121	3	1	1	67	20	1	0.693147181
121	7	2	1	67	20	1	0.693147181
121	7	3	1	67	20	1	0.693147181
121	7	4	1	67	20	1	0.693147181
122	41	0	1	41	22	0	2.079441542
122	4	1	1	41	22	1	0.693147181
122	18	2	1	41	22	1	0.693147181
122	2	3	1	41	22	1	0.693147181
122	5	4	1	41	22	1	0.693147181
123	12	0	0	12	31	0	2.079441542
123	6	1	0	12	31	1	0.693147181
123	4	2	0	12	31	1	0.693147181
123	0	3	0	12	31	1	0.693147181
123	2	4	0	12	31	1	0.693147181
124	7	0	1	7	28	0	2.079441542
124	2	1	1	7	28	1	0.693147181
124	1	2	1	7	28	1	0.693147181
124	1	3	1	7	28	1	0.693147181
124	0	4	1	7	28	1	0.693147181
126	52	0	0	52	42	0	2.079441542
126	40	1	0	52	42	1	0.693147181
126	20	2	0	52	42	1	0.693147181
126	23	3	0	52	42	1	0.693147181
126	12	4	0	52	42	1	0.693147181
128	22	0	1	22	23	0	2.079441542
128	0	1	1	22	23	1	0.693147181
128	2	2	1	22	23	1	0.693147181
128	4	3	1	22	23	1	0.693147181
128	0	4	1	22	23	1	0.693147181
129	13	0	1	13	40	0	2.079441542
129	5	1	1	13	40	1	0.693147181
129	4	2	1	13	40	1	0.693147181
129	0	3	1	13	40	1	0.693147181
129	3	4	1	13	40	1	0.693147181
130	23	0	0	23	37	0	2.079441542
130	5	1	0	23	37	1	0.693147181
130	6	2	0	23	37	1	0.693147181
130	6	3	0	23	37	1	0.693147181
130	5	4	0	23	37	1	0.693147181
135	10	0	0	10	28	0	2.079441542
135	14	1	0	10	28	1	0.693147181
135	13	2	0	10	28	1	0.693147181
135	6	3	0	10	28	1	0.693147181
135	0	4	0	10	28	1	0.693147181
137	46	0	1	46	33	0	2.079441542
137	11	1	1	46	33	1	0.693147181
137	14	2	1	46	33	1	0.693147181
137	25	3	1	46	33	1	0.693147181
137	15	4	1	46	33	1	0.693147181
139	36	0	1	36	21	0	2.079441542
139	10	1	1	36	21	1	0.693147181
139	5	2	1	36	21	1	0.693147181
139	3	3	1	36	21	1	0.693147181
139	8	4	1	36	21	1	0.693147181
141	52	0	0	52	36	0	2.079441542
141	26	1	0	52	36	1	0.693147181
141	12	2	0	52	36	1	0.693147181
141	6	3	0	52	36	1	0.693147181
141	22	4	0	52	36	1	0.693147181
143	38	0	1	38	35	0	2.079441542
143	19	1	1	38	35	1	0.693147181
143	7	2	1	38	35	1	0.693147181
143	6	3	1	38	35	1	0.693147181
143	7	4	1	38	35	1	0.693147181
145	33	0	0	33	24	0	2.079441542
145	12	1	0	33	24	1	0.693147181
145	6	2	0	33	24	1	0.693147181
145	8	3	0	33	24	1	0.693147181
145	4	4	0	33	24	1	0.693147181
147	7	0	1	7	25	0	2.079441542
147	1	1	1	7	25	1	0.693147181
147	1	2	1	7	25	1	0.693147181
147	2	3	1	7	25	1	0.693147181
147	3	4	1	7	25	1	0.693147181
201	18	0	0	18	23	0	2.079441542
201	4	1	0	18	23	1	0.693147181
201	4	2	0	18	23	1	0.693147181
201	6	3	0	18	23	1	0.693147181
201	2	4	0	18	23	1	0.693147181
202	42	0	0	42	36	0	2.079441542
202	7	1	0	42	36	1	0.693147181
202	9	2	0	42	36	1	0.693147181
202	12	3	0	42	36	1	0.693147181
202	14	4	0	42	36	1	0.693147181
203	36	0	1	36	26	0	2.079441542
203	6	1	1	36	26	1	0.693147181
203	10	2	1	36	26	1	0.693147181
203	8	3	1	36	26	1	0.693147181
203	8	4	1	36	26	1	0.693147181
204	11	0	1	11	25	0	2.079441542
204	2	1	1	11	25	1	0.693147181
204	1	2	1	11	25	1	0.693147181
204	0	3	1	11	25	1	0.693147181
204	0	4	1	11	25	1	0.693147181
205	87	0	0	87	26	0	2.079441542
205	16	1	0	87	26	1	0.693147181
205	24	2	0	87	26	1	0.693147181
205	10	3	0	87	26	1	0.693147181
205	9	4	0	87	26	1	0.693147181
206	50	0	0	50	26	0	2.079441542
206	11	1	0	50	26	1	0.693147181
206	0	2	0	50	26	1	0.693147181
206	0	3	0	50	26	1	0.693147181
206	5	4	0	50	26	1	0.693147181
208	22	0	1	22	32	0	2.079441542
208	4	1	1	22	32	1	0.693147181
208	3	2	1	22	32	1	0.693147181
208	2	3	1	22	32	1	0.693147181
208	4	4	1	22	32	1	0.693147181
209	41	0	1	41	25	0	2.079441542
209	8	1	1	41	25	1	0.693147181
209	6	2	1	41	25	1	0.693147181
209	5	3	1	41	25	1	0.693147181
209	7	4	1	41	25	1	0.693147181
210	18	0	0	18	28	0	2.079441542
210	0	1	0	18	28	1	0.693147181
210	0	2	0	18	28	1	0.693147181
210	3	3	0	18	28	1	0.693147181
210	3	4	0	18	28	1	0.693147181
211	32	0	1	32	35	0	2.079441542
211	1	1	1	32	35	1	0.693147181
211	3	2	1	32	35	1	0.693147181
211	1	3	1	32	35	1	0.693147181
211	5	4	1	32	35	1	0.693147181
213	111	0	0	111	31	0	2.079441542
213	37	1	0	111	31	1	0.693147181
213	29	2	0	111	31	1	0.693147181
213	28	3	0	111	31	1	0.693147181
213	29	4	0	111	31	1	0.693147181
214	56	0	1	56	21	0	2.079441542
214	18	1	1	56	21	1	0.693147181
214	11	2	1	56	21	1	0.693147181
214	28	3	1	56	21	1	0.693147181
214	13	4	1	56	21	1	0.693147181
215	18	0	0	18	32	0	2.079441542
215	3	1	0	18	32	1	0.693147181
215	5	2	0	18	32	1	0.693147181
215	2	3	0	18	32	1	0.693147181
215	5	4	0	18	32	1	0.693147181
217	20	0	0	20	21	0	2.079441542
217	3	1	0	20	21	1	0.693147181
217	0	2	0	20	21	1	0.693147181
217	6	3	0	20	21	1	0.693147181
217	7	4	0	20	21	1	0.693147181
218	24	0	1	24	41	0	2.079441542
218	6	1	1	24	41	1	0.693147181
218	3	2	1	24	41	1	0.693147181
218	4	3	1	24	41	1	0.693147181
218	0	4	1	24	41	1	0.693147181
219	12	0	0	12	29	0	2.079441542
219	3	1	0	12	29	1	0.693147181
219	4	2	0	12	29	1	0.693147181
219	3	3	0	12	29	1	0.693147181
219	4	4	0	12	29	1	0.693147181
220	9	0	0	9	21	0	2.079441542
220	3	1	0	9	21	1	0.693147181
220	4	2	0	9	21	1	0.693147181
220	3	3	0	9	21	1	0.693147181
220	4	4	0	9	21	1	0.693147181
221	16	0	1	16	32	0	2.079441542
221	3	1	1	16	32	1	0.693147181
221	5	2	1	16	32	1	0.693147181
221	4	3	1	16	32	1	0.693147181
221	3	4	1	16	32	1	0.693147181
222	17	0	0	17	32	0	2.079441542
222	2	1	0	17	32	1	0.693147181
222	3	2	0	17	32	1	0.693147181
222	3	3	0	17	32	1	0.693147181
222	5	4	0	17	32	1	0.693147181
225	22	0	1	22	26	0	2.079441542
225	1	1	1	22	26	1	0.693147181
225	23	2	1	22	26	1	0.693147181
225	19	3	1	22	26	1	0.693147181
225	8	4	1	22	26	1	0.693147181
226	28	0	0	28	25	0	2.079441542
226	8	1	0	28	25	1	0.693147181
226	12	2	0	28	25	1	0.693147181
226	2	3	0	28	25	1	0.693147181
226	8	4	0	28	25	1	0.693147181
227	55	0	0	55	30	0	2.079441542
227	18	1	0	55	30	1	0.693147181
227	24	2	0	55	30	1	0.693147181
227	76	3	0	55	30	1	0.693147181
227	25	4	0	55	30	1	0.693147181
228	25	0	1	25	21	0	2.079441542
228	2	1	1	25	21	1	0.693147181
228	3	2	1	25	21	1	0.693147181
228	0	3	1	25	21	1	0.693147181
228	1	4	1	25	21	1	0.693147181
230	9	0	0	9	40	0	2.079441542
230	2	1	0	9	40	1	0.693147181
230	1	2	0	9	40	1	0.693147181
230	2	3	0	9	40	1	0.693147181
230	1	4	0	9	40	1	0.693147181
232	13	0	1	13	36	0	2.079441542
232	0	1	1	13	36	1	0.693147181
232	0	2	1	13	36	1	0.693147181
232	0	3	1	13	36	1	0.693147181
232	0	4	1	13	36	1	0.693147181
234	10	0	0	10	19	0	2.079441542
234	3	1	0	10	19	1	0.693147181
234	1	2	0	10	19	1	0.693147181
234	4	3	0	10	19	1	0.693147181
234	2	4	0	10	19	1	0.693147181
236	12	0	1	12	37	0	2.079441542
236	1	1	1	12	37	1	0.693147181
236	4	2	1	12	37	1	0.693147181
236	3	3	1	12	37	1	0.693147181
236	2	4	1	12	37	1	0.693147181
238	47	0	0	47	22	0	2.079441542
238	13	1	0	47	22	1	0.693147181
238	15	2	0	47	22	1	0.693147181
238	13	3	0	47	22	1	0.693147181
238	12	4	0	47	22	1	0.693147181




<<code>>
seizure<- read.csv('C:/Users/UX303LB/Desktop/class/R/seizure.csv')
head(seizure)
glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)

#第一題
y<-(seizure$y)
X<-cbind(1,seizure$trt,seizure$age)
head(X)
time<-exp(seizure$ltime)

f3<-function(beta){
  mu<-exp(X%*%beta)*time
  ymu<-y-mu
  gradient<-t(X)%*%ymu
  hessian<--t(X)%*%diag(as.vector(mu))%*%X
  return(list(gradient, hessian))
}

newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
  x <- x0
  f3.x <- f3(x)
  n <- 0
  while ((max(abs(f3.x[[1]])) > tol) & (n < n.max)) {
    x <- x - solve.default(f3.x[[2]], tol= 1e-19)%*% f3.x[[1]]
    f3.x <- f3(x)
    n <- n + 1


  }
  if (n == n.max) {
    cat('newton failed to converge\n')
  } else {
    return(x)
  }
}

newton(f3,c(0, 0, 0),1e-9)

#第二題
result<-glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)
b<-newton(f3,c(0, 0, 0),1e-9)
vcov(result)
solve(-hessian)
solve(-f3(b)[[2]], tol= 1e-19)


#第三題
beta<-newton(f3,c(0, 0, 0),1e-9)
result<-glm(seizure$y~seizure$trt+seizure$age,offset=seizure$ltime,family=poisson)
logLik(result)
loglike<-sum((-exp(X%*%beta)%*%time)+(y%*%X%*%beta)+(y%*%(seizure$ltime))-(log(factorial(y))))

a<--exp(X%*%beta)%*%time
b<-y%*%X%*%beta
c<-y%*%(seizure$ltime)
d<--log(factorial(y))
loglike<-sum(a+diag(b)+diag(c)+d)

loglike<-sum(-mu+diag(y%*%log(mu))-log(factorial(y)))
loglike<-sum(-mu+diag(y%*%log(mu))-log(factorial(y)))