Definimos las funciones a usar.
Valor_F <- function (Probabilidad,Numerador,Denominador){
Resultado <- qf(Probabilidad,Numerador,Denominador,lower.tail = FALSE)
return (Resultado)
}
Valor_t <- function (Probabilidad, Libertad){
Resultado <- qt(Probabilidad,Libertad,lower.tail = FALSE)
return (Resultado)
}
Valor_Absoluto <- function (Numero){
if (Numero<0){
return (Numero*(-1))
}else{
return (Numero)
}
}
M_X <- function (Vector_X1,Vector_X2,Vector_X3,Vector_X4,Vector_X5,Vector_X6,Vector_X7,Vector_X8,Vector_X9,Vector_X10,Vector_X11,Vector_X12,Vector_X13,Vector_X14,Vector_X15){
Vector <- c()
for (i in 1:length(Vector_X1)){
Vector <- c(Vector,1)
}
M <- cbind(Vector,Vector_X1,Vector_X2,Vector_X3,Vector_X4,Vector_X5,Vector_X6,Vector_X7,Vector_X8,Vector_X9,Vector_X10,Vector_X11,Vector_X12,Vector_X13,Vector_X14,Vector_X15)
return (M)
}
H <- function (M_X){
T_M_X <- t(M_X)
A1 <- T_M_X%*%M_X
A <- solve(A1)
return (M_X%*%A%*%T_M_X)
}
Matriz_b <- function(M_X,M_Y){
T_M_X <- t(M_X)
A1 <- T_M_X%*%M_X
A <- solve(A1)
B <- T_M_X%*%M_Y
Resultado <- A%*%B
row.names(Resultado) <- c("b0","b1","b2","b3","b4","b5","b6","b7","b8","b9","b10","b11","b12","b13","b14","b15")
colnames(Resultado) <- c("Coeficientes de Regresion")
return (Resultado)
}
Matriz_Valores_Ajustados <- function (M_X,M_Y){
Resultado <- H(M_X)%*%M_Y
colnames(Resultado) <- c("Valores ajustados")
return (Resultado)
}
Matriz_residuos <- function (M_X,M_Y){
Matriz_H <- H(M_X)
Matriz_Identidad <- diag(1, nrow = nrow(Matriz_H), ncol = ncol(Matriz_H))
A <- Matriz_Identidad-Matriz_H
Resultado <- A%*%M_Y
colnames(Resultado) <- c("Residuos")
return (Resultado)
}
Matriz_Xh <- function (X_h1,X_h2,X_h3,X_h4,X_h5,X_h6,X_h7,X_h8,X_h9,X_h10,X_h11,X_h12,X_h13,X_h14,X_h15){
Vector <- c(1,X_h1,X_h2,X_h3,X_h4,X_h5,X_h6,X_h7,X_h8,X_h9,X_h10,X_h11,X_h12,X_h13,X_h14,X_h15)
Vector2 <- cbind(Vector)
colnames(Vector2) <- c("Valores Xh")
return (Vector2)
}
Matriz_Gorro_Yh <- function (M_X,M_Y,M_Xh){
b <- Matriz_b(M_X,M_Y)
a <- t(M_Xh)
resultado <- a%*%b
colnames(resultado) <- c("Respuesta media estimada")
return(resultado)
}
SSTO <- function (M_Y){
n <- length(M_Y)
Matriz_Identidad <- diag(1, nrow = n, ncol = n)
Matriz_1 <- matrix(1,nrow = n, ncol = n)
Transpuesta <- t(M_Y)
A <- (Matriz_Identidad-(1/n)*Matriz_1)
return (Transpuesta%*%A%*%M_Y)
}
SSE <- function (M_X,M_Y){
n <- length(M_Y)
Matriz_H <- H(M_X)
Matriz_Identidad <- diag(1, nrow = nrow(Matriz_H), ncol = ncol(Matriz_H))
Transpuesta <- t(M_Y)
A <- Matriz_Identidad - Matriz_H
return (Transpuesta%*%A%*%M_Y)
}
SSR <- function (M_X,M_Y){
n <- length(M_Y)
Matriz_H <- H(M_X)
Matriz_1 <- matrix(1, nrow = n, ncol= n)
Transpuesta <- t(M_Y)
A <- (Matriz_H - (1/n)*Matriz_1)
return (Transpuesta%*%A%*%M_Y)
}
ANOVA <- function (M_X,M_Y,Probabilidad){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
ssr <- SSR(M_X,M_Y)
sse <- SSE(M_X,M_Y)
ssto <- SSTO(M_Y)
msr <- ssr/p_1
mse <- sse/n_p
Estadistico_Prueba <- msr/mse
Critico <- Valor_F(Probabilidad,p_1,n_p)
FUENTE_VARIACION <- c("Regresion","Error","Total")
SS <- c(ssr,sse,ssto)
GL <- c(p_1,n_p,n-1)
MS <- c(msr,mse,"-")
ESTADISTICO_PRUEBA <- c(Estadistico_Prueba,"-","-")
VALOR_CRITICO <- c(Critico,"-","-")
M <- cbind(FUENTE_VARIACION,SS,GL,MS,ESTADISTICO_PRUEBA,VALOR_CRITICO)
return (M)
}
Matriz_Varianza_Covarianza_de_Matriz_Residuo <- function (M_X,M_Y){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
mse <- SSE(M_X,M_Y)/n_p
Matriz_H <- H(M_X)
Matriz_Identidad <- diag(1, nrow = nrow(Matriz_H), ncol = ncol(Matriz_H))
A <- Matriz_Identidad-Matriz_H
return (mse[1]*A)
}
COEFICIENTE_DETERMINACION <- function (M_X,M_Y){
ssr <- SSR(M_X,M_Y)
ssto <- SSTO(M_Y)
resultado <- ssr/ssto
colnames(resultado) <- c("Coeficiente de determinacion")
row.names(resultado) <- c("[1]")
return (resultado)
}
Contribucion <- function (Matriz_X_Original,Matriz_X_Nueva,y){
E_Full <- ANOVA(Matriz_X_Original,y,0.05)[2,2]
E_Red <- ANOVA(Matriz_X_Nueva,y,0.05)[2,2]
return (as.numeric(E_Full) - as.numeric(E_Red))
}
#El numero de columnas y filas dependen del numero de variables, si tenemos 2 variables entonces tenemos 3 columnas y 3 filas
S_b <- function (M_X,M_Y){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado <- (mse[1]*A)^(1/2)
colnames(resultado) <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)
row.names(resultado) <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)
return (resultado)
}
S_Gorro_Yh <- function (M_X,M_Y,M_Xh){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado_a <- mse[1]*A
Transpuesta_M_Xh <- t(M_Xh)
return ((Transpuesta_M_Xh%*%resultado_a%*%M_Xh)^(1/2))
}
#EL VALOR DE PROBAILIDAD INGRESAR ES ALPHA/2, PARAMETROS TIENE EL VALOR
#DEL NUMERO DE PARAMETROS A CALCULAR LIMITE
Valor_Bonferroni <- function (M_X,M_Y,Probabilidad,Parametros){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
t <- Valor_t(Probabilidad*(1/Parametros),n_p)
}
#INTERVALO DE CONFIANZA CONJUNTO DE BONFERRONI
#EL VALOR DEL ESTIMADO i ES EL B_i CORRESPONDIENTE, PARA B_1 ES 1
#EL VALOR DE PROBABILIDAD ES ALPHA/2
Intervalo_Confianza_Conjunta_Bonferroni <- function (M_X,M_Y,Probabilidad,Parametros,Estimado_i){
B <- Valor_Bonferroni(M_X,M_Y,Probabilidad,Parametros)
b_i <- Matriz_b(M_X,M_Y)[Estimado_i+1]
s_b <- S_b(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Izquierda <- b_i - B*s_b
Derecha <- b_i + B*s_b
Resultado <- c(Izquierda,"<Betta<",Derecha)
return (Resultado)
}
Intervalo_Confianza_E_Yh <- function (M_X,M_Y,M_Xh,Probabilidad){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
Gorro_Yh <- Matriz_Gorro_Yh(M_X,M_Y,M_Xh)[1]
t <- Valor_t(Probabilidad,n_p)
s_gorro_yh <- S_Gorro_Yh(M_X,M_Y,M_Xh)[1]
Izquierda <- Gorro_Yh - t*s_gorro_yh
Derecha <- Gorro_Yh + t*s_gorro_yh
Resultado <- c(Izquierda,"<E{Yh}<",Derecha)
return (Resultado)
}
Prueba_Hipotesis_Bk <- function (M_X,M_Y,Probabilidad,Estimado_i){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
M1 <- "Rechazamos Ho"
M2 <- "No rechazamos Ho"
Critico <- Valor_t(Probabilidad,n_p)
b_i <- Matriz_b(M_X,M_Y)[Estimado_i+1]
s_b <- S_b(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Estadistico_Prueba <- b_i/s_b
print("Estadistico Prueba")
print(Valor_Absoluto(Estadistico_Prueba))
print("Critico")
print(Critico)
if (Valor_Absoluto(Estadistico_Prueba)>Critico){
return (M1)
}else{
return (M2)
}
}
Prueba_Hipotesis_Inferencia_Variable<-function (Modelo_Original_X,Modelo_Nuevo_X,M_Y,Probabilidad){
M1 <- "Las variables describen bien el experimento"
M2 <- "Las variables no describen bien el experimento"
sse_modelo_original <- as.numeric(ANOVA(Modelo_Original_X,M_Y,Probabilidad)[2,2])
sse_modelo_nuevo <- as.numeric(ANOVA(Modelo_Nuevo_X,M_Y,Probabilidad)[2,2])
gl_modelo_original <- as.numeric(ANOVA(Modelo_Original_X,M_Y,Probabilidad)[2,3])
gl_modelo_nuevo <- as.numeric(ANOVA(Modelo_Nuevo_X,M_Y,Probabilidad)[2,3])
A <- (sse_modelo_original-sse_modelo_nuevo)/(gl_modelo_original-gl_modelo_nuevo)
B <- (sse_modelo_nuevo/gl_modelo_nuevo)
Estadistico_Prueba <- A/B
Critico <- Valor_F(Probabilidad,gl_modelo_original-gl_modelo_nuevo,gl_modelo_nuevo)
print("Estadistico de Prueba")
print(Estadistico_Prueba)
print("Valor Critico")
print(Critico)
if (Estadistico_Prueba>Critico){
return(M1)
}else{
return(M2)
}
}
Datos del problema y valor de Alpha
#VALOR DE ALPHA
Alpha <- 0.05
#DATOS DEL PROBLEMA
datos <- read.table("ApendiceC.02.txt")
data <- cbind(datos[,4:9],datos[,11:16])
categoricos <- datos[,17]
cat <- matrix(0, nrow=length(categoricos), ncol=3)
for(i in 1:length(categoricos)){
if(categoricos[i]==1){
cat[i,1] <- 1
}else if(categoricos[i]==2){
cat[i,2] <- 1
}else if(categoricos[i]==3){
cat[i,3] <- 1
}
}
data <- cbind(data,cat)
data <- as.matrix(data)
colnames(data) <- c("Land_area","Total_population","%population_aged_18-34","%population_aged_65+",
"Active_physicians","Hospital_beds","%highschool_graduates",
"%bachelor_degrees","%below_poverty_level","%unemployment","Per_capita_income",
"Total_personal_income","NE_region","NC_region","S_region")
x1 <- data[,1]
x2 <- data[,2]
x3 <- data[,3]
x4 <- data[,4]
x5 <- data[,5]
x6 <- data[,6]
x7 <- data[,7]
x8 <- data[,8]
x9 <- data[,9]
x10 <- data[,10]
x11 <- data[,11]
x12 <- data[,12]
x13 <- data[,13]
x14 <- data[,14]
x15 <- data[,15]
Matriz_Y <- cbind(datos[,10])
Matriz_X <- M_X(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15)
M_Xh <- Matriz_Xh(1150,1346597,45.1,12.9,10560,4680,78.9,34.1,8.4,13.5,42653,11468,0,1,0)
Analisis del modelo completo
#MATRIZ DE B
Matriz_b(Matriz_X,Matriz_Y)
## Coeficientes de Regresion
## b0 -4.464259e+04
## b1 -3.155656e+00
## b2 2.335731e-01
## b3 2.238672e+02
## b4 1.032252e+01
## b5 -5.381499e+00
## b6 3.602665e+00
## b7 -2.459025e+02
## b8 1.146088e+02
## b9 1.161689e+03
## b10 -2.154378e+02
## b11 2.322607e+00
## b12 -7.042099e+00
## b13 -4.402755e+03
## b14 -4.571025e+03
## b15 -2.347399e+03
#CALCULAR EL ANOVA
ANOVA(Matriz_X,Matriz_Y,Alpha)
## FUENTE_VARIACION SS GL MS
## [1,] "Regresion" "1256460785665.33" "15" "83764052377.6885"
## [2,] "Error" "232454753160.525" "424" "548242342.359728"
## [3,] "Total" "1488915538825.85" "439" "-"
## ESTADISTICO_PRUEBA VALOR_CRITICO
## [1,] "152.786543296079" "1.69002155308825"
## [2,] "-" "-"
## [3,] "-" "-"
#PRUEBAD DE HIPOTESIS PARA Bk
for (L1 in 0:15){
print(paste0("La prueba de Hipotesis para B",L1," es "))
Prueba_Hipotesis_Bk(Matriz_X,Matriz_Y,Alpha/2,L1)
}
## [1] "La prueba de Hipotesis para B0 es "
## [1] "Estadistico Prueba"
## [1] 1.199719
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B1 es "
## [1] "Estadistico Prueba"
## [1] 3.479707
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B2 es "
## [1] "Estadistico Prueba"
## [1] 11.46377
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B3 es "
## [1] "Estadistico Prueba"
## [1] 0.5252393
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B4 es "
## [1] "Estadistico Prueba"
## [1] 0.02669341
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B5 es "
## [1] "Estadistico Prueba"
## [1] 1.682203
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B6 es "
## [1] "Estadistico Prueba"
## [1] 1.554221
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B7 es "
## [1] "Estadistico Prueba"
## [1] 0.6962332
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B8 es "
## [1] "Estadistico Prueba"
## [1] 0.2957849
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B9 es "
## [1] "Estadistico Prueba"
## [1] 2.367579
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B10 es "
## [1] "Estadistico Prueba"
## [1] 0.3147388
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B11 es "
## [1] "Estadistico Prueba"
## [1] 3.306594
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B12 es "
## [1] "Estadistico Prueba"
## [1] 7.417388
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B13 es "
## [1] "Estadistico Prueba"
## [1] 0.9959947
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B14 es "
## [1] "Estadistico Prueba"
## [1] 1.086315
## [1] "Critico"
## [1] 1.965575
## [1] "La prueba de Hipotesis para B15 es "
## [1] "Estadistico Prueba"
## [1] 0.5670516
## [1] "Critico"
## [1] 1.965575
#CALCULAR EL INTERVALO DE CONFIANZA CONJUNTA BONFERRONI
#EL VALOR DE Parametros ES DEL NUMERO DE PARAMETROS A ESTIMAR AL MISMO TIEMPO.
Parametros = 15
for (L2 in 0:15){
print(paste0("El intervalo de confianza Bonferroni para B",L2," es "))
print(Intervalo_Confianza_Conjunta_Bonferroni(Matriz_X,Matriz_Y,Alpha/2,Parametros,L2))
}
## [1] "El intervalo de confianza Bonferroni para B0 es "
## [1] "-154486.374431241" "<Betta<" "65201.2035274886"
## [1] "El intervalo de confianza Bonferroni para B1 es "
## [1] "-5.8326829174948" "<Betta<" "-0.478628843278868"
## [1] "El intervalo de confianza Bonferroni para B2 es "
## [1] "0.17342788600882" "<Betta<" "0.293718341180023"
## [1] "El intervalo de confianza Bonferroni para B3 es "
## [1] "-1034.30177227487" "<Betta<" "1482.03612166148"
## [1] "El intervalo de confianza Bonferroni para B4 es "
## [1] "-1131.20775710667" "<Betta<" "1151.85279444836"
## [1] "El intervalo de confianza Bonferroni para B5 es "
## [1] "-14.8249500055566" "<Betta<" "4.06195195213338"
## [1] "El intervalo de confianza Bonferroni para B6 es "
## [1] "-3.23986698999669" "<Betta<" "10.445197516668"
## [1] "El intervalo de confianza Bonferroni para B7 es "
## [1] "-1288.49356675822" "<Betta<" "796.688532569511"
## [1] "El intervalo de confianza Bonferroni para B8 es "
## [1] "-1029.18489751966" "<Betta<" "1258.40252100679"
## [1] "El intervalo de confianza Bonferroni para B9 es "
## [1] "-286.720073562336" "<Betta<" "2610.09839666489"
## [1] "El intervalo de confianza Bonferroni para B10 es "
## [1] "-2236.0234972083" "<Betta<" "1805.14792802339"
## [1] "El intervalo de confianza Bonferroni para B11 es "
## [1] "0.249123129353592" "<Betta<" "4.39609035164024"
## [1] "El intervalo de confianza Bonferroni para B12 es "
## [1] "-9.84467131097153" "<Betta<" "-4.23952644598319"
## [1] "El intervalo de confianza Bonferroni para B13 es "
## [1] "-17451.6352190619" "<Betta<" "8646.12498321689"
## [1] "El intervalo de confianza Bonferroni para B14 es "
## [1] "-16992.222481583" "<Betta<" "7850.17264837508"
## [1] "El intervalo de confianza Bonferroni para B15 es "
## [1] "-14567.366991824" "<Betta<" "9872.56906897423"
#CALCULAR EL INTERVALO DE CONFIANZA PARA E{Yh}
Intervalo_Confianza_E_Yh(Matriz_X,Matriz_Y,M_Xh,Alpha/2)
## [1] "160683.30010021" "<E{Yh}<" "322539.178513033"
#MATRIZ DE COVARIANZA DEL ERROR
# DEJAMOS EN COMENTARIO ESTA FUNCION PORQUE ARROJA MAS DE 100 ELEMENTOS.
#Matriz_Varianza_Covarianza_de_Matriz_Residuo(Matriz_X,Matriz_Y)
Analisis modelo 1
#FUNCION
M_X_3 <- function (Vector_X1,Vector_X2,Vector_X3){
Vector <- c()
for (i in 1:length(Vector_X1)){
Vector <- c(Vector,1)
}
M <- cbind(Vector,Vector_X1,Vector_X2,Vector_X3)
return (M)
}
#APARTADO PARA ANALIZAR EL MODELO 1
MX_C <- M_X_3(x13,x14,x15)
#COEFICIENTE DE DETERMINACION
COEFICIENTE_DETERMINACION(MX_C,Matriz_Y)
## Coeficiente de determinacion
## [1] 0.01205893
#MATRIZ DE B
Matriz_b_C <- function(M_X,M_Y){
T_M_X <- t(M_X)
A1 <- T_M_X%*%M_X
A <- solve(A1)
B <- T_M_X%*%M_Y
Resultado <- A%*%B
row.names(Resultado) <- c("b0","b13","b14","b15")
colnames(Resultado) <- c("Coeficientes de Regresion")
return (Resultado)
}
Matriz_b_C(MX_C,Matriz_Y)
## Coeficientes de Regresion
## b0 40211.73
## b13 -17125.70
## b14 -18430.99
## b15 -13220.75
#ANOVA
ANOVA(MX_C,Matriz_Y,Alpha)
## FUENTE_VARIACION SS GL MS
## [1,] "Regresion" "17954726745.0336" "3" "5984908915.01122"
## [2,] "Error" "1470960812080.82" "436" "3373763330.46059"
## [3,] "Total" "1488915538825.85" "439" "-"
## ESTADISTICO_PRUEBA VALOR_CRITICO
## [1,] "1.77395635934965" "2.62536382495787"
## [2,] "-" "-"
## [3,] "-" "-"
#PRUEBA DE HIPOTESIS PARA BK
S_b_C <- function (M_X,M_Y){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado <- (mse[1]*A)^(1/2)
colnames(resultado) <- c(1,2,3,4)
row.names(resultado) <- c(1,2,3,4)
return (resultado)
}
Prueba_Hipotesis_Bk_C <- function (M_X,M_Y,Probabilidad,Estimado_i){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
M1 <- "Rechazamos Ho"
M2 <- "No rechazamos Ho"
Critico <- Valor_t(Probabilidad,n_p)
b_i <- Matriz_b_C(M_X,M_Y)[Estimado_i+1]
s_b <- S_b_C(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Estadistico_Prueba <- b_i/s_b
print("Estadistico Prueba")
print(Valor_Absoluto(Estadistico_Prueba))
print("Critico")
print(Critico)
if (Valor_Absoluto(Estadistico_Prueba)>Critico){
return (M1)
}else{
return (M2)
}
}
#B0
Prueba_Hipotesis_Bk_C(MX_C,Matriz_Y,Alpha/2,0)
## [1] "Estadistico Prueba"
## [1] 6.074923
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#B13
Prueba_Hipotesis_Bk_C(MX_C,Matriz_Y,Alpha/2,1)
## [1] "Estadistico Prueba"
## [1] 1.957126
## [1] "Critico"
## [1] 1.96542
## [1] "No rechazamos Ho"
#B14
Prueba_Hipotesis_Bk_C(MX_C,Matriz_Y,Alpha/2,2)
## [1] "Estadistico Prueba"
## [1] 2.127466
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#B15
Prueba_Hipotesis_Bk_C(MX_C,Matriz_Y,Alpha/2,3)
## [1] "Estadistico Prueba"
## [1] 1.627228
## [1] "Critico"
## [1] 1.96542
## [1] "No rechazamos Ho"
#INTERVALO BONFERRONI
Intervalo_Confianza_Conjunta_Bonferroni_C <- function (M_X,M_Y,Probabilidad,Parametros,Estimado_i){
B <- Valor_Bonferroni(M_X,M_Y,Probabilidad,Parametros)
b_i <- Matriz_b_C(M_X,M_Y)[Estimado_i+1]
s_b <- S_b_C(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Izquierda <- b_i - B*s_b
Derecha <- b_i + B*s_b
Resultado <- c(Izquierda,"<Betta<",Derecha)
return (Resultado)
}
#B0
Intervalo_Confianza_Conjunta_Bonferroni_C(MX_C,Matriz_Y,Alpha/2,3,0)
## [1] "24303.8747152764" "<Betta<" "56119.5798301781"
#B13
Intervalo_Confianza_Conjunta_Bonferroni_C(MX_C,Matriz_Y,Alpha/2,3,1)
## [1] "-38155.2102659571" "<Betta<" "3903.8139729297"
#B14
Intervalo_Confianza_Conjunta_Bonferroni_C(MX_C,Matriz_Y,Alpha/2,3,2)
## [1] "-39251.2175236684" "<Betta<" "2389.24445969536"
#B15
Intervalo_Confianza_Conjunta_Bonferroni_C(MX_C,Matriz_Y,Alpha/2,3,3)
## [1] "-32746.4936874721" "<Betta<" "6304.98651043866"
#VALORES PROPUESTOS
Matriz_Xh_C <- function (X_h1,X_h2,X_h3){
Vector <- c(1,X_h1,X_h2,X_h3)
Vector2 <- cbind(Vector)
colnames(Vector2) <- c("Valores Xh")
return (Vector2)
}
Matriz_Gorro_Yh_C <- function (M_X,M_Y,M_Xh){
b <- Matriz_b_C(M_X,M_Y)
a <- t(M_Xh)
resultado <- a%*%b
colnames(resultado) <- c("Respuesta media estimada")
return(resultado)
}
S_Gorro_Yh_C <- function (M_X,M_Y,M_Xh){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado_a <- mse[1]*A
Transpuesta_M_Xh <- t(M_Xh)
return ((Transpuesta_M_Xh%*%resultado_a%*%M_Xh)^(1/2))
}
Intervalo_Confianza_E_Yh_C <- function (M_X,M_Y,M_Xh,Probabilidad){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
Gorro_Yh <- Matriz_Gorro_Yh_C(M_X,M_Y,M_Xh)[1]
t <- Valor_t(Probabilidad,n_p)
s_gorro_yh <- S_Gorro_Yh_C(M_X,M_Y,M_Xh)[1]
Izquierda <- Gorro_Yh - t*s_gorro_yh
Derecha <- Gorro_Yh + t*s_gorro_yh
Resultado <- c(Izquierda,"<E{Yh}<",Derecha)
return (Resultado)
}
M_Xh_C <- Matriz_Xh_C(1,0,0)
Intervalo_Confianza_E_Yh_C(MX_C,Matriz_Y,M_Xh_C,Alpha/2)
## [1] "11837.5441430344" "<E{Yh}<" "34334.5141093928"
Analisis modelo 2
MX_J <- M_X_3(x3,x10,x12)
#COEFICIENTE DE DETERMINACION
COEFICIENTE_DETERMINACION(MX_J,Matriz_Y)
## Coeficiente de determinacion
## [1] 0.7187298
#MATRIZ DE B
Matriz_b_J <- function(M_X,M_Y){
T_M_X <- t(M_X)
A1 <- T_M_X%*%M_X
A <- solve(A1)
B <- T_M_X%*%M_Y
Resultado <- A%*%B
row.names(Resultado) <- c("b0","b3","b10","b12")
colnames(Resultado) <- c("Coeficientes de Regresion")
return (Resultado)
}
Matriz_b_J(MX_J,Matriz_Y)
## Coeficientes de Regresion
## b0 -38730.680144
## b3 754.527639
## b10 2172.418884
## b12 3.806716
#ANOVA
ANOVA(MX_J,Matriz_Y,Alpha)
## FUENTE_VARIACION SS GL MS
## [1,] "Regresion" "1070127919005.58" "3" "356709306335.192"
## [2,] "Error" "418787619820.279" "436" "960522063.80798"
## [3,] "Total" "1488915538825.85" "439" "-"
## ESTADISTICO_PRUEBA VALOR_CRITICO
## [1,] "371.370236849138" "2.62536382495787"
## [2,] "-" "-"
## [3,] "-" "-"
#PRUEBA DE HIPOTESIS PARA BK
S_b_J <- function (M_X,M_Y){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado <- (mse[1]*A)^(1/2)
colnames(resultado) <- c(1,2,3,4)
row.names(resultado) <- c(1,2,3,4)
return (resultado)
}
Prueba_Hipotesis_Bk_J <- function (M_X,M_Y,Probabilidad,Estimado_i){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
M1 <- "Rechazamos Ho"
M2 <- "No rechazamos Ho"
Critico <- Valor_t(Probabilidad,n_p)
b_i <- Matriz_b_J(M_X,M_Y)[Estimado_i+1]
s_b <- S_b_J(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Estadistico_Prueba <- b_i/s_b
print("Estadistico Prueba")
print(Valor_Absoluto(Estadistico_Prueba))
print("Critico")
print(Critico)
if (Valor_Absoluto(Estadistico_Prueba)>Critico){
return (M1)
}else{
return (M2)
}
}
#B0
Prueba_Hipotesis_Bk_J(MX_J,Matriz_Y,Alpha/2,0)
## [1] "Estadistico Prueba"
## [1] 3.093224
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#B3
Prueba_Hipotesis_Bk_J(MX_J,Matriz_Y,Alpha/2,1)
## [1] "Estadistico Prueba"
## [1] 2.049015
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#B10
Prueba_Hipotesis_Bk_J(MX_J,Matriz_Y,Alpha/2,2)
## [1] "Estadistico Prueba"
## [1] 3.297393
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#B12
Prueba_Hipotesis_Bk_J(MX_J,Matriz_Y,Alpha/2,3)
## [1] "Estadistico Prueba"
## [1] 33.07055
## [1] "Critico"
## [1] 1.96542
## [1] "Rechazamos Ho"
#INTERVALO BONFERRONI
Intervalo_Confianza_Conjunta_Bonferroni_J <- function (M_X,M_Y,Probabilidad,Parametros,Estimado_i){
B <- Valor_Bonferroni(M_X,M_Y,Probabilidad,Parametros)
b_i <- Matriz_b_J(M_X,M_Y)[Estimado_i+1]
s_b <- S_b_J(M_X,M_Y)[Estimado_i+1,Estimado_i+1]
Izquierda <- b_i - B*s_b
Derecha <- b_i + B*s_b
Resultado <- c(Izquierda,"<Betta<",Derecha)
return (Resultado)
}
#B0
Intervalo_Confianza_Conjunta_Bonferroni_J(MX_J,Matriz_Y,Alpha/2,3,0)
## [1] "-68822.1518695482" "<Betta<" "-8639.2084186831"
#B3
Intervalo_Confianza_Conjunta_Bonferroni_J(MX_J,Matriz_Y,Alpha/2,3,1)
## [1] "-130.44461196668" "<Betta<" "1639.49989029422"
#B10
Intervalo_Confianza_Conjunta_Bonferroni_J(MX_J,Matriz_Y,Alpha/2,3,2)
## [1] "589.085306582347" "<Betta<" "3755.75246176841"
#B12
Intervalo_Confianza_Conjunta_Bonferroni_J(MX_J,Matriz_Y,Alpha/2,3,3)
## [1] "3.53007990621227" "<Betta<" "4.08335183926461"
#VALORES PROPUESTOS
Matriz_Xh_J <- function (X_h1,X_h2,X_h3){
Vector <- c(1,X_h1,X_h2,X_h3)
Vector2 <- cbind(Vector)
colnames(Vector2) <- c("Valores Xh")
return (Vector2)
}
Matriz_Gorro_Yh_J <- function (M_X,M_Y,M_Xh){
b <- Matriz_b_J(M_X,M_Y)
a <- t(M_Xh)
resultado <- a%*%b
colnames(resultado) <- c("Respuesta media estimada")
return(resultado)
}
S_Gorro_Yh_J <- function (M_X,M_Y,M_Xh){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
sse <- SSE(M_X,M_Y)
mse <- sse/n_p
Transpuesta <- t(M_X)
A <- solve(Transpuesta%*%M_X)
resultado_a <- mse[1]*A
Transpuesta_M_Xh <- t(M_Xh)
return ((Transpuesta_M_Xh%*%resultado_a%*%M_Xh)^(1/2))
}
Intervalo_Confianza_E_Yh_J <- function (M_X,M_Y,M_Xh,Probabilidad){
n <- length(M_Y)
p_1 <- ncol(M_X)-1
n_p <- (n-1)-p_1
Gorro_Yh <- Matriz_Gorro_Yh_J(M_X,M_Y,M_Xh)[1]
t <- Valor_t(Probabilidad,n_p)
s_gorro_yh <- S_Gorro_Yh_J(M_X,M_Y,M_Xh)[1]
Izquierda <- Gorro_Yh - t*s_gorro_yh
Derecha <- Gorro_Yh + t*s_gorro_yh
Resultado <- c(Izquierda,"<E{Yh}<",Derecha)
return (Resultado)
}
M_Xh_J <- Matriz_Xh_J(62.4,32.65,79462)
Intervalo_Confianza_E_Yh_J(MX_J,Matriz_Y,M_Xh_J,Alpha/2)
## [1] "332465.006947639" "<E{Yh}<" "431076.148627504"