# data[].attributes.body
halfScore = sum(pred.prob >= 0.5)

confM <- confusionMatrix(data=as.factor(pred.val),reference=as.factor(model_data$PoolsLocal))

freqtab = data.frame(Predicted=table(pred.val),Observed = table(model_data$PoolsLocal))
freqtab = select(freqtab,-Observed.Var1)

# data[].attributes.body
# Exercise 2c: Make a ROC curve (and learn how to write a loop).

# As we change the cut value τ
# we can observe a trade-off between the true positive rate (TPR)
# and the false positive rate (FPR).
# A higher cut value will result in fewer positive predictions (both true and false),
# and therefore a lower TPR and FPR.
# On the other hand, a lower cut value will result
# in more positive predictions (both true and false),
# and therefore a higher TPR and FPR.
# This is what the ROC curve illustrates.

# Calculate the TPR and the FPR for three cut values: 0, 0.35 and 1.
# Here, I do it for the first cut value by relying on indexation (Hermansen 2023, ch 3).

#Empty objects
TPR=NULL;FPR=NULL

#Cut values
tau<-c(0,0.35,1)

#Indexation
df$preds0<-predict(mod1,df,"response") > tau[1]

truePos<-sum(df$PoolsLocal==1&df$preds0==1,na.rm = T)
Pos<-sum(df$PoolsLocal,na.rm = T)

falsePos<-sum(df$PoolsLocal==0&df$preds0==1,na.rm = T)
Neg<-sum(df$PoolsLocal==0,na.rm = T)

#Indexation
TPR[1]<-truePos/Pos
FPR[1]<-falsePos/Neg

# You can do it for the following two cut values by
# replacing 1 by 2, and then 3 in ([1]).

#Indexation
df$preds0<-predict(mod1,df,"response") > tau[2]

truePos<-sum(df$PoolsLocal==1&df$preds0==1,na.rm = T)
Pos<-sum(df$PoolsLocal,na.rm = T)

falsePos<-sum(df$PoolsLocal==0&df$preds0==1,na.rm = T)
Neg<-sum(df$PoolsLocal==0,na.rm = T)

#Indexation
TPR[2]<-truePos/Pos
FPR[2]<-falsePos/Neg

#Indexation
df$preds0<-predict(mod1,df,"response") > tau[3]

truePos<-sum(df$PoolsLocal==1&df$preds0==1,na.rm = T)
Pos<-sum(df$PoolsLocal,na.rm = T)

falsePos<-sum(df$PoolsLocal==0&df$preds0==1,na.rm = T)
Neg<-sum(df$PoolsLocal==0,na.rm = T)

#Indexation
TPR[3]<-truePos/Pos
FPR[3]<-falsePos/Neg

# Show values
TPR
FPR

# Make a plot where you define FPR along the x-axis
# and the TPR along the y-axis.
# Draw a line between the three dots

#Create plot with labeled axes
plot(FPR, TPR,
type="l",
xlab="False Positive Rate (FPR)",
ylab="True Positive Rate (TPR)",
main = "ROC plot for mod1")

# Use the R code from the first ROC curve,
# wrap it in a for loop, and let the tau vary from 0 to 1 through increments of 0.05.
# Use the code to illustrate the ROC curve of the second model in this problem sheet

# This code uses a for loop to iterate over each tau value in the tau vector.
# Inside the loop, the code calculates the true positive rate (TPR)
# and false positive rate (FPR) for the current tau value
# and saves them in the TPR and FPR vectors, respectively.
# The loop continues until all tau values have been processed.

#Empty objects
TPR=NULL;FPR=NULL

#Cut values
tau_seq <- seq(from = 0, to = 1, by = 0.05)

for (i in 1:length(tau_seq)) {

# Indexation
df$preds0 <- predict(mod2, df, "response") > tau_seq[i]

truePos <- sum(df$PoolsLocal == 1 & df$preds0 == 1, na.rm = T)
Pos <- sum(df$PoolsLocal, na.rm = T)

falsePos <- sum(df$PoolsLocal == 0 & df$preds0 == 1, na.rm = T)
Neg <- sum(df$PoolsLocal == 0, na.rm = T)

# Indexation
TPR[i] <- truePos / Pos
FPR[i] <- falsePos / Neg
}

TPR
FPR

#Create plot with labeled axes
plot(FPR, TPR,
type="l",
xlab="False Positive Rate (FPR)",
ylab="True Positive Rate (TPR)",
main = "ROC plot for mod2")

# data[].attributes.body
# Exercise 2b: True positives vs. false positives

# Create a cross table where you compare
# your predicted 0s and 1s with the actual observed values.
# This is called a confusion matrix.

# Rename variables
observed <- df$PoolsLocal
predicted_mod1 <- tau_cut_numeric1
predicted_mod2 <- tau_cut_numeric2

# create table
conf_mat1 <- table(observed, predicted_mod1)
conf_mat1

conf_mat2 <- table(observed, predicted_mod2)
conf_mat2

# Calculate the true positive rate (sensitivity),
# defined as the number of correctly predicted 1s divided by all observed 1s

# convert confusion matrix to proportions, where each row sums to 1.
conf_mat1 %>%
prop.table(., margin = 1)
# sensitivity is seen the in lower right corner (68,8%)

conf_mat2 %>%
prop.table(., margin = 1)
# sensitivity is seen the in lower right corner (76%)

# Calculate the false positive rate,
# defined as the number of incorrectly predicted 0s divided by the number of observed 0s

# Mod1: the false positive rate is seen in the upper right corner (45,5%)
# Mod2: the false positive rate is seen in the upper right corner (37%)

# This means that with the selected tau cut of 0.35, mod2 perfoms best

#Alternative
#Cut value;tau
tau <- 0.35

#Predicted  outcome
df$preds1  <- predict(mod1, df, "response")>tau
df$preds2  <- predict(mod2, df, "response")>tau

#True positives
truePos1 <- sum(df$PoolsLocal == 1 & df$preds1 == 1, na.rm = T)
truePos2 <- sum(df$PoolsLocal == 1 & df$preds2 == 1, na.rm = T)

#All positives
Pos<-sum(df$PoolsLocal,na.rm = T)

#The rate
truePos1/Pos
truePos2/Pos

# data[].attributes.body
# Exercise 2a: Predicted values

# Calculate the predicted values from model 1 and 2.
pred_values1 <- predict(mod1,
#In-sample prediction
newdata = df,
#I want probabilities
type = "response")

pred_values2 <- predict(mod2,
#In-sample prediction
newdata = df,
#I want probabilities
type = "response")

# How many observations have a predicted probability above 0.5?
# count number of observations with predicted probability above 0.5,
# excluding missing values
sum(pred_values1 > 0.5 & complete.cases(pred_values1)) # 0 observations
sum(pred_values2 > 0.5 & complete.cases(pred_values2)) # 172 observations

# Decide on a τ, the cut point beyond which
# the predicted probability translates to a 1. Why did you choose that one?

# I use the base-line model (intercept only model / proportion of 1s).
# It is the value at which I’d predict as much wrong
# in one direction as the other without any other predictors in the model.

# τ = 0.35
# (i.e the proportion of 1s in the dataset is 35%)

# However, the correct threshold depends on the cost of false negatives and false positives

# Create a frequency table.
# How many 0s and 1s did you predict?
# How does it compare with the observed frequency of PoolsLocal?

# create logical vector based on tau cut of 0.35
tau_cut1 <- pred_values1 >= 0.35
tau_cut2 <- pred_values2 >= 0.35

# convert logical vector to numeric vector with 1s and 0s
tau_cut_numeric1 <- as.numeric(tau_cut1)
tau_cut_numeric2 <- as.numeric(tau_cut2)

# create frequency table
table(tau_cut_numeric1)
table(tau_cut_numeric2)

# Observed frequency of Poolslocal
table(df$PoolsLocal)

# Both model 1 and model 2 predict a higher number of 1s (pools ressources)
# than the observed outcome.

# data[].attributes.body
#ROC(K OUT) curves

#Empty objects
TPR=NULL;
FPR=NULL

#Cut values
tau<-c(0,0.35,1)

#Indexation
for(i in 1:length(tau)) {
preds0<-predict(mod1,df,"response") > tau[i]
truePos<-sum(df$PoolsLocal==1 & preds0==1,na.rm = T)
Pos<-sum(df$PoolsLocal,na.rm = T)

falsePos<-sum(df$PoolsLocal==0 & preds0==1,na.rm = T)
Neg<-sum(df$PoolsLocal==0,na.rm = T)

#Indexation
TPR[i]<-truePos/Pos
FPR[i]<-falsePos/Neg

}
roc <- data.frame(TPR, FPR)
roc %>%
ggplot() +
geom_point(aes(x=FPR, y=TPR)) +
geom_segment(aes(x=FPR[1], xend=FPR[2], y=TPR[1], yend=TPR[2])) +
geom_segment(aes(x=FPR[2], xend=FPR[3], y=TPR[2], yend=TPR[3]))

# data[].attributes.body
TPR1 = NULL
FPR1 = NULL
tau_loop <- seq(from=0, to=1, by=0.05)

for(j in 1:length(tau_loop)) {

preds0<-predict(mod1,df,"response") > tau_loop[j]
truePos<-sum(df$PoolsLocal==1 & preds0==1,na.rm = T)
Pos<-sum(df$PoolsLocal,na.rm = T)

falsePos<-sum(df$PoolsLocal==0 & preds0==1,na.rm = T)
Neg<-sum(df$PoolsLocal==0,na.rm = T)

#Indexation
TPR1[j] <- truePos / Pos
FPR1[j] <- falsePos / Neg
}

# data[].attributes.body
#count obs - model1=0, model2=172
tau<-0.5
predmatrix %>%
filter(predmod1>tau) %>%
count()

# I have no idea how to choose a proper cutoff, so I just worked with  .3.

#generate vars
tau <- 0.3
predmatrix<-predmatrix %>%
mutate(preds1 = if_else(predmod1>tau, 1, 0),
preds2 = if_else(predmod2>tau, 1, 0))

#confusion matrix as prop table
prop.table(table(predmatrix$df.PoolsLocal, predmatrix$preds1))
prop.table(table(predmatrix$df.PoolsLocal, predmatrix$preds2))
#model 2 slightly better

#TPR = TP / TP+FN
#FPR = FP / FP+TN
tpr1 <- 172/(172+78)
fpr1 <- 208/(208+249)
tpr2 <- 210/(210+36)
fpr2 <- 207/(207+233)

# data[].attributes.body
###### Exercise 2a

# How many observations have a predicted probability above 0.5?

sum(df$preds1 > 0.5)
# 0 values

sum(df$preds2 > 0.5, na.rm = T)
# 172 values.

# Decide on a τ, the cut point beyond which the predicted probability translates to a 1.
# Why did you choose that one?

# I choose the mean of df$PoolsLocal.  I choose this cut-off because it represents the proportion of 0 and 1 in df$PoolsLocal.
I looked up other ways to decide on a cut-off value and found Youden J-value, but I was not able to calculate a J-value. Hence, i went with the mean of df$PoolsLocal.

# Now, I find the mean of the values
summary(df$PoolsLocal)

# Recode preds as 1 or 0 depending on whether predicted value is above of below the mean
# i.e. Tao = mean of PoolsLocal.
df$pred1 = ifelse(df$preds1 > mean(df$PoolsLocal),
yes = 1,
no = 0)
df$pred2 = ifelse(df$preds2 > mean(df$PoolsLocal),
yes = 1,
no = 0)

# Create a confusion matrix with the observed and predicted value.
df[, c("PoolsLocal", "pred1")] %>%
table %>%
prop.table
# 59% of observation are correctly classified as 0 and 1.

df[, c("PoolsLocal", "pred2")] %>%
table %>%
prop.table
# 67% of observations are correctly classified as 0 and 1.

###### Exercise 2b

# Normal frequency table to calculate True Positive Rate and False Positive Rate:
table(df$PoolsLocal, df$pred1)
TRP1 = 172 / (172+78)
FRP1 = 208 / (208+249)

table(df$PoolsLocal, df$pred2)
TRP2 = 187 / (187+59)
FRP2 = 163 / (163+277)