Расчет вероятности покрытия интервала прогнозирования (PICP)

 В цифровом картографировании почв (DSM) мы делаем прогнозы пространственного распределения свойств почвы, которое сопряжено с неопределенностями / ошибками. Для количественной оценки точности мы разбиваем данные на обучающий и тестовый наборы, где мы обучаем модель машинного обучения (например, случайный лес, аддитивные модели, сплайны и т. д.) на обучающем наборе и оцениваем модель на тестовом наборе. Затем мы предсказываем ковариаты, чтобы составить непрерывную карту, а также предсказываем неопределенности (например, стандартное отклонение, стандартная ошибка, предельный диапазон предсказания и т. д.). Но как мы оцениваем оценки неопределенности, которые мы только что предсказали? Одним из вариантов является расчет вероятности покрытия интервала прогнозирования (PICP), который измеряет, насколько хорошо прогнозы попадают в сформулированный интервал прогнозирования (PI). PI похожи на доверительные интервалы, но используют прогнозируемые данные для определения интервалов.

Эта функция будет в пакете rafikisol R, однако я еще не загрузил ее, так как есть еще несколько функций, которые я хочу добавить. Тем не менее, мы превратим его в функцию под названием calcPICP(), которая была в основном взята из книги «Использование R для цифрового картирования почвы» Malone et al., (2017).

Функция принимает 3 параметра:  x = a data frame of the data, response = the vector of the measured data (e.g., data$response), and pred = the predicted value (e.g., data$predicted).

calcPICP = function(data, response, pred){

#We сначала получаем остатки модели

res = response - pred

#Then получаем стандартное отклонение остатков и объединяем с данными.

data$stdev = sd(res)

#We чем сделать ряд квантилей из нормального кумулятивного распределения.

qp <- qnorm(c(0.995, 0.9875, 0.975, 0.95, 0.9, 0.8, 0.7, 0.6, 0.55, 0.525))

#Then составляем матрицу с длиной строки данных и столбцов qp

vMat <- matrix(NA, nrow = nrow(data), ncol = length(qp))

#Now мы должны обойти квантили и умножить его на стандартное отклонение, чтобы получить серию стандартных ошибок с разными интервалами предсказания.



}
#Make другую матрицу, такую же, как и раньше, для верхних пределов 

uMat <- matrix(NA, nrow = nrow(data), ncol = length(qp))

#We вычислить верхние пределы, добавив ряд стандартных ошибок к предсказаниям модели.

uMat[,  i] <- pred + vMat[, i]
}
#We составляем еще одну матрицу для нижних пределов #We вычисляем нижние пределы


, вычитая ряд из предсказанных значений.

lMat[, i] <- pred - vMat[, i] }
#Now мы хотим увидеть, какие интервалы прогнозирования охватывают измеренные данные, создавая матрицу из 1 и 0.

bMat[,  i]<-as.numeric(response <= uMat[,  i]  &
response >= lMat[, i])
}

#To вычисляем PICP, берем colsums/nrow*100 для матрицы 0s и 1s


#Make вектор уровней достоверности 

#We помещенный в кадр данных для построения графика 

#Since мы хотим, чтобы PICP к CI был линией 1:1, мы также вычисляем коэффициент корреляции соответствия Лина (CCC) с пакетом критериев R.


#Make имя правильно


#must добавить значения осей для построения

ccc$y = 90 #y axis

мы можем построить график PICP в CI, добавить линию 1:1 и CCC geom_point()+ 

#add точки geom_text(data = ccc,aes(x= x, y =y, label = paste("CCC = ",round(CCC, 2))))+ #add значение CCC

geom_point()+ #add points
geom_text(data = ccc,aes(x= x, y =y, label = paste("CCC = ",round(CCC, 2))))+ #add CCC value
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = 'red')+ #add 1:1 line
labs(x = ‘Confidence level’, y = “PICP”, title = “PICP to confidence level”)+ #labels
theme_bw() #make it look good

#make это выглядит хорошо

, #Now мы хотим вернуть список графика, а также кадр данных общих результатов.

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

vMat[, i] <- data$stdev * qp[i]

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

lMat <- matrix(NA, nrow = nrow(data), ncol = length(qp))

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

bMat <- matrix(NA, nrow = nrow(data), ncol = length(qp))

for(i in 1:ncol(bMat)){

picp <- colSums(bMat)/nrow(bMat)*100

cl <- c(99, 97.5, 95, 90, 80, 60, 40, 20, 10, 5)

results <- data.frame(picp = picp, cl = cl)

ccc <- as.data.frame(yardstick::ccc_vec(results$picp, results$cl))

names(ccc) = "CCC" #name

ccc$x = 10 #x axis

p = ggplot(data = results, aes(x= cl, y = picp))+ #add data

return(setNames(list(p, results), c("Plot", "Results")))

}

Теперь у нас есть функция, дающая нам график преобразования PICP в CI и результаты. Это полезно при запуске многих моделей, и теперь мы можем просто подключить данные.



Теперь постройте данные.