#craga de datos y estimacion del modelo

options(scipen = 999999)
library(wooldridge)
library(dplyr)
library(stargazer)
data("lawsch85")
lawsch85<-na.omit(lawsch85)
modelo_salario<-lm(formula =lsalary~LSAT+GPA+libvol+cost+rank,data = lawsch85)
stargazer(modelo_salario,title = 'modelo de regresion salario',type = "html")
modelo de regresion salario
Dependent variable:
lsalary
LSAT 0.007
(0.005)
GPA 0.221**
(0.111)
libvol 0.0002**
(0.0001)
cost 0.00001
(0.00000)
rank -0.003***
(0.0004)
Constant 8.871***
(0.680)
Observations 90
R2 0.837
Adjusted R2 0.827
Residual Std. Error 0.111 (df = 84)
F Statistic 85.997*** (df = 5; 84)
Note: p<0.1; p<0.05; p<0.01

Definición funciones para el cálculo de Theil y su descomposición, debe estar instalada la librería DescTools:

#Bias Proportion
Um<-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado) 
}
#Variance Proportion
Us<-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc<-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U<-function(pronosticado,observado){
   library(DescTools)
  RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}

scrib

options(scipen = 999999) 
library(dplyr) 
library(caret) 
library(DescTools) 
library(stargazer) 
set.seed(100) 
numero_de_muestras<-175
muestras<- lawsch85$lsalary %>%
  createDataPartition(p = 0.75,
                      times = numero_de_muestras,
                      list = TRUE)
# Listas vacias, que contendran los datos de entrenamiento, los pronosticos para los datos de prueba, y para las estadisticas de cada muestra
Modelos_Entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance<-vector(mode = "list",
                              length = numero_de_muestras)
#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- lawsch85[muestras[[j]], ]
Datos_Prueba<- lawsch85[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data=Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]]<-data.frame( 
            R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
                    Datos_Entrenamiento$lsalary),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$lsalary),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$lsalary),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$lsalary)*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$lsalary,type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$lsalary),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$lsalary),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$lsalary)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$lsalary),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$lsalary),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$lsalary),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$lsalary)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$lsalary,
                         type = 1), # También se puede usar la función que creamos: THEIL_U
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$lsalary),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$lsalary),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$lsalary)
            )
}

resultados

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo",
            type = "html",
            digits = 3)
Medidas de Performance Datos del Modelo
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 175 0.835 0.014 0.808 0.825 0.844 0.877
RMSE 175 0.108 0.005 0.091 0.105 0.111 0.117
MAE 175 0.086 0.004 0.074 0.084 0.089 0.094
MAPE 175 0.815 0.036 0.701 0.791 0.840 0.890
THEIL 175 0.005 0.0002 0.004 0.005 0.005 0.006
Um 175 0.000 0.000 0 0 0 0
Us 175 0.046 0.004 0.033 0.043 0.049 0.054
Uc 175 0.969 0.004 0.960 0.966 0.972 0.981 
bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación",
            type = "html",
            digits = 3)
Medidas de Performance Simulación
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 175 0.824 0.053 0.633 0.795 0.860 0.946
RMSE 175 0.113 0.016 0.072 0.101 0.124 0.160
MAE 175 0.091 0.014 0.062 0.081 0.101 0.138
MAPE 175 0.865 0.130 0.591 0.765 0.952 1.294
THEIL 175 0.005 0.001 0.003 0.005 0.006 0.008
Um 175 0.054 0.068 0.00000 0.006 0.088 0.418
Us 175 0.082 0.089 0.00000 0.013 0.118 0.528
Uc 175 0.914 0.113 0.468 0.868 1.000 1.052