-
Notifications
You must be signed in to change notification settings - Fork 0
/
logistic_regression_and_gof_functions.r
105 lines (100 loc) · 3.53 KB
/
logistic_regression_and_gof_functions.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
99
100
101
102
103
104
105
######
# makes the sucess array, that is if there are 5 dead out 10 at 0.1 mg it creates
# 1 1 1 1 1 0 0 0 0 0
# 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
# parameter "df" is a dataframe object
######
MakeBinomialPairs <- function(df){
retX <- list()
retY <- list()
for(i in 1:nrow(df)){
x <- rep(log(df$dose[i]), df$batch[i])
y <- c(rep(1, df$dead[i]), rep(0, df$batch[i]-df$dead[i]))
retX <- append(retX, x)
retY <- append(retY, y)
retX <- unlist(retX, recursive=TRUE)
retY <- unlist(retY, recursive=TRUE)
}
return(myDf <- data.frame(logDose=retX, dead=retY))
}
####
# A function to do the Hosmer-Lemeshow test in R.
# Original by www.r-bloggers.com/author/denishaine/
# modified by [email protected]
# param y are the observed and yhat the adjusted by the model and g?
####
Hosmerlem <- function (y, yhat, g)
{
try
#todo what to do if cut is not unique maybe add a try loop.
cutyhat <- tryCatch({
cut(yhat, breaks = quantile(yhat, probs = seq(0, 1, 1/g)), include.lowest = T)
},
error=function(cond) {
message("(╯°□°）╯ ︵ ┻━┻ The model has many ties in its predicted probabilities")
message("(ಠ_ಠ) Hosmer-Lemeshow test is not applicable in this case")
message("(ಠ_ಠ) Too few covariate values?")
message("Here's the original error message:")
message(cond)
# Choose a return value in case of error
stop()
},
warning=function(cond) {
message("Here's the original warning message:")
message(cond)
# Choose a return value in case of warning
},
finally={
}
)
obs <- xtabs(cbind(1 - y, y) ~ cutyhat)
expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
chisq <- sum((obs - expect)^2/expect)
P <- 1 - pchisq(chisq, g - 2)
c("X^2" = chisq, Df = g - 2, "P(>Chi)" = P)
}
######
#reverse prediction of x given a logistic regression
# parameter "y" is the observed value, intercept is the beta 0 and the sloep is beta 1
######
PredictX <- function(y, intercept, slope){
x <- (log(y/(1-y)) - intercept) / slope
return(x)
}
#####################################################################
# Summary: calculate the CI for beta1 "slope/dose" parameter using the delta method
# based on code by Dason Kurkiewicz http://dasonk.com/r/2013/02/09/Using-the-delta-method/
# modified by [email protected] Nov 25th, 2014
# parameters
################################################
ConfInt <- function(estimatedDose, p, glmModel){
standerr <- deltamethod(~ (log(p/(1-p)) - x1)/x2, coef(glmModel), vcov(glmModel))
ci <- estimatedDose + c(-1, 1) * 1.96 * standerr
return(ci)
}
######
# Checks if a library is present. If present it loads.
# If the library is not present tries to install
# and load the requested library.
# [email protected] 25th Nov, 2014
# parameter(s) "libname" is the name of the library to check in string
######
Check_n_Load <- function(libName){
if(require(libName, character.only=TRUE)){
print(paste("@(^_^)@ ", libName," is loaded correctly"))
}
else {
print(paste("@(◔ ◡ ◔)@, trying to install", libName ))
install.packages(libName)
if(require(libName)){
print(paste(libName, "@(^_^)@ installed and loaded"))
} else {
stop(paste("(╯°□°）╯ ︵ ┻━┻ could not install", libname))
}
}
}
Check_n_Load("msm")# needed for deltamethod function
Check_n_Load("rms")# needed le Cessie – van Houwelingen – Copas – Hosmer unweighted sum of squares test
# for global goodness of fit