-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMAST30027Assignment1.R
103 lines (71 loc) · 2.76 KB
/
MAST30027Assignment1.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
load("assign1.Robj")
model = glm(pima_subset$test~pima_subset$bmi,family = binomial(logit))
summary(model)
##install.packages('package_name', dependencies=TRUE, repos='http://cran.rstudio.com/')
#Q1a
library(faraway)
data(orings)
logL<- function(beta, orings){
eta = cbind(rep(1,length(orings$temp)),orings$temp)%*%beta
return(sum(orings$damage*log(pnorm(eta))+(6-orings$damage)*log(1-pnorm(eta))))
}
#maximize
(betahat = optim(c(10,-.1), logL,orings = orings, control = list(fnscale=-1))$par)
##Q1b
beta0 = betahat[1]
beta1 = betahat[2]
(eta = beta0+orings$temp*beta1)
#cdf of standard normal
D0 = pnorm(eta)
#first derivative(pdf of standard normal)
D1 = dnorm(eta)
#second derivative
D2=1/sqrt(2*pi)*exp(-(eta^2)/2)*(-1*eta)
x = orings$temp
y = orings$damage
m = 6
(Lb1b1=sum( (y*x^2*(D0*D2-D1*D1))/(D0^2)+(y-m)*(x^2*((1-D0)*D2+D1*D1))/((1-D0)*(1-D0))))
(Lb1b0=sum( y*x*(D0*D2-D1*D1)/(D0^2)+(y-m)*x*((1-D0)*D2+D1*D1)/((1-D0)*(1-D0))))
(Lb0b0=sum( y*(D0*D2-D1*D1)/(D0^2)+(y-m)*((1-D0)*D2+D1*D1)/((1-D0)*(1-D0))))
(J = -matrix(c(Lb0b0,Lb1b0,Lb1b0,Lb1b1),2,2))
(variance = solve(J))
(Beta0CI = c(betahat[1] - qnorm(0.975)*sqrt(variance[1,1]),betahat[1] + qnorm(0.975)*sqrt(variance[1,1])))
(Beta1CI = c(betahat[2] - qnorm(0.975)*sqrt(variance[2,2]),betahat[2] + qnorm(0.975)*sqrt(variance[2,2])))
##Q1c Likelihood Ratio Test
#phatN for the null model
n=rep(6,length(y))
(phatN = sum(y)/sum(n))
#Max log likelihood for Reduced model
(MaxlogL.R = sum(y*log(phatN))+sum(6-y)*log(1-phatN))
#Max log likelihood for Full model
(MaxlogL.F = logL(betahat, orings))
#Likelihood Ratio Statistics
(LR = -2*(MaxlogL.R-MaxlogL.F))
#LR follow chi-square distribution with 1 degree of freedom
(pValue = pchisq(df=1,LR,lower.tail = FALSE))
#Since p-value is less than 0.05,reject null hypothesis (H0)
#Thus beta-1 is not equal to 0
#Q1d
t=c(1,31)
#estimate
eta = t(t)%*%betahat
(pnorm(eta))
#the estimate of damage when temperature equals 31 Fahrenheit is 0.9896084
#95%CI for estimate
etaCI = c(t%*%betahat-qnorm(0.975)*sqrt(t(t)%*%variance%*%t),
t%*%betahat+qnorm(0.975)*sqrt(t(t)%*%variance%*%t))
(pnorm(etaCI))
#The 95% Confidence interval for estimate probability of damage when temperature
#equals to 31 Fahrenheit is (0.7169598, 0.9999744)
#Q1e
#Plot the data
plot(damage/6 ~ temp, orings, xlim=c(25,85), ylim=c(0,1),xlab="Temperature", ylab="Prob of damage")
x <- seq(25,85,1)
#use the glm with logit link to fit the model
logitmod <- glm(cbind(damage,6-damage) ~ temp, family=binomial, orings)
(betahatLogit=logitmod$coefficients)
lines(x, ilogit(betahatLogit[1] + betahatLogit[2]*x), col="red",)
#use probit link to fit the model
lines(x, pnorm(betahat[1] + betahat[2]*x), col="green")
#add legend
legend(60,1,legend = c("logit model","probit model"), col = c("red","green"),lty=1)