How to do a Calibration plot by cohorts

Question

I have been trying to do a Calibration plot by the variable cohorts, instead of doing a single one for all of my sample in order to see if some cohort underestimate and others overestimate. However, I can not manage to do it and I keep getting errors. Could you please help me how to make it work or perhaps there is a different library in R that allows me to do it.

Thank you so much

My code is the following:

<code>
data = SCORE_DSAExtended

#Codigo vara borrar los casos con missing values
data_pp = na.omit(data)

# Create cvd10 column based on conditions
data_pp$cvd10 <- ifelse(data_pp$cvd_rip == 0 & data_pp$survtime > 10, 0,
                        ifelse(data_pp$cvd_rip == 1 & data_pp$survtime <= 10, 1,
                               ifelse(data_pp$cvd_rip == 1 & data_pp$survtime > 10, 0, NA)))

# Discard cases where cvd_rip = 0 and survtime < 10
data_pp <- data_pp[complete.cases(data_pp$cvd10), ]


#Create variable blood pressure (BP)
library(dplyr)

data_pp <- data_pp %>%
  mutate(bp = case_when(
    sbp < 120 & dbp < 80 ~ "Normal",
    sbp >= 120 & sbp <= 139 | dbp >= 80 & dbp <= 89 ~ "Elevada",
    sbp >= 140 & sbp <= 159 | dbp >= 90 & dbp <= 99 ~ "Grado1",
    TRUE ~ "Grado2"
  ))

bp = ordered(data_pp$bp, levels=c("Normal", "Elevada", "Grado1", "Grado2"))

data_pp$bp <- factor(data_pp$bp, levels = c("Normal", "Elevada", "Grado1", "Grado2"))

#Fitting logistic regression

#Variables as factors
data_pp$smoker = as.factor(data_pp$smoker)
data_pp$cohort = as.factor(data_pp$cohort)
data_pp$sex = as.factor(data_pp$sex)
data_pp$gl_int = as.factor(data_pp$gl_int)
data_pp$prev_mi = as.factor(data_pp$prev_mi)

lreg2 = glm(cvd10~I(age/10)+cohort+smoker+sex+gl_int+prev_mi+hdliu+choliu+I(choliu^2)+bp+sex:smoker+sex:gl_int+I(age/10):smoker, data=data_pp, family=binomial(link="logit"))

pred2 =  predict(lreg2, newdata = data_pp, type = "response")

#Calibration Belt as Measure of Calibration
library(givitiR)
cb = givitiCalibrationBelt(data_pp$cvd10, pred2, devel="internal")
plot(cb)
<\code>

In this last bit how do I include that I want it by cohort. I have also tried some other functions available in R, that are not working, like:

<code>
calibration_plot(data = data_pp, obs = "cvd10", pred = "pred2", y_lim = c(0, 1), x_lim=c(0, 1),title = "Calibration plot", group = "cohort")

cohorts=data_pp$cohort
calib = calibrate(pred2, cohorts, method="boot", B=200)
plot(calib, main="Calibration Plot by Cohort")
<\code>