# install.packages("Car")
library(car)
## Cargando paquete requerido: carData
#install.packages("Cairo")
library(Cairo)

#install.packages("readxl")
library(readxl)


#install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)
## Cargando paquete requerido: xts
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Adjuntando el paquete: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
#install.packages("ggplot2")
library(ggplot2)

#install.packages("lmtest")
library(lmtest)
#install.packages("dplyr")
library(dplyr)
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ✔ readr     2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first()  masks xts::first()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::last()   masks xts::last()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ISLR)

Base de datos

La base de datos que se trabaja es de biomasa del tallo de las platas de café, relacionando las variables altura del tallo frenta a peso del tallo

datos<-read_excel("C:/Users/Acer/Downloads/taller.xlsx")
datos
## # A tibble: 102 × 5
##    Grosor_inferior Grosor_superior Altura  Peso clase_de_cafe
##              <dbl>           <dbl>  <dbl> <dbl> <chr>        
##  1            12.9             5.9   204.   854 bourbon      
##  2            11.5             5.9   177.   670 castillo     
##  3            11.8             5.9   181.   722 castillo     
##  4            11.5             4.5   206.   817 castillo     
##  5            13.8             4.8   240.   653 bourbon      
##  6            13.3             5.9   168.   848 castillo     
##  7            13               6.3   176.   939 castillo     
##  8            10.4             3     182.   619 castillo     
##  9            12.6             5.1   175.   790 castillo     
## 10            10.1             6.7   188.   702 castillo     
## # ℹ 92 more rows

Gráfica de dispersión

plot(datos$Altura, datos$Peso, main = "Dispersion entre Altura del tallo y Peso del tallo", xlab = "Altura del tallo", ylab = "Peso del tallo")

Se puede observar que las variables altura del tallo y peso del tallo tiene una relación directa, ya que a medida que aumenta la altura también lo hace el peso

Correlación

cor_test <- cor.test(datos$Altura, datos$Peso)
cor_test
## 
##  Pearson's product-moment correlation
## 
## data:  datos$Altura and datos$Peso
## t = 7.9818, df = 100, p-value = 2.488e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.488630 0.729766
## sample estimates:
##       cor 
## 0.6238251

Las variables altura y peso tiene una correlacion apriximada de 0,62, esta siendo mayor a 0,5 se considera que las variables tiene una correlacion positiva fuerte.

Modelo de Regresión

regresion <- lm(Peso ~ Altura , data = datos)
summary(regresion)
## 
## Call:
## lm(formula = Peso ~ Altura, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -568.95 -218.73  -35.42  204.85  843.28 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -200.0400   154.2900  -1.297    0.198    
## Altura         5.9126     0.7408   7.982 2.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 290.6 on 100 degrees of freedom
## Multiple R-squared:  0.3892, Adjusted R-squared:  0.383 
## F-statistic: 63.71 on 1 and 100 DF,  p-value: 2.488e-12

\[y=\beta_0 + \beta_1(x)+e \]

Con el análisis de regresión encontramos los parámetros\(\beta_0\) y \(\beta_1\), donde el modelo queda de la siguiente manera, tiene una variación en la altura del tallo en un 43,5%

\[y= -200.04 + 5.91(x) \]

A continuación se logra observar la predicción de modelo frente a los datos reales de las medidas de los tallos de los arboles de café.

#library(ggplot2)
ggplot(datos, aes(x=Altura, y=Peso)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth(method=lm)
## `geom_smooth()` using formula = 'y ~ x'

En promedio se considera que cuando aumenta la altura de la planta tambien lo hace el peso del tallo 7 gramos, sin embargo, no se puede relacionar con la realidad un tallo que tenga una altura de cero(0) y a su vez un peso de -449.91. Dado que la regresion estima como logramos ver, estimaciones mas alla de los valores observados.

plot(datos$Altura, datos$Peso, xlab='Altura del tallo', ylab='Peso del tallo')
abline(regresion, col = "red", lwd = 2)

# Independecia de errores. A contiación se muestra graficamente el comportamiento de los residuos frente a los valores ajustados, donde se observa que los valores se distribuye de manera aleatoria al rededor de cero(0) sin formar patrones evidentes.

plot(regresion$fitted.values, regresion$residuals,
     main = 'Residuos vs Ajustados',
     xlab = 'Valores Ajustados',
     ylab = 'Residuos',
     pch = 19)
abline(h = 0, col = "red", lwd = 2)

De acuerdo a lo anterior es necesario aplicar el test estadístico de Durbin-Watson que permite dectectar la autocorrelación de los residuos.

dwtest(regresion)
## 
##  Durbin-Watson test
## 
## data:  regresion
## DW = 2.3293, p-value = 0.953
## alternative hypothesis: true autocorrelation is greater than 0

De acuerdo al estadístico de Durbin-Watson el valor obtenido se aproximó a 2 y el p-valor es mayor a 0.05 lo que indica que no se puede rechazar la hipotesis nula (H0) de independencia de los residuos, por lo tanto, se cumple el supuesto de independencia.

Homocedasticidad

Se aplica la prueba de Breusch-Pagan para verificar si la varianza de los errores es constante.

bptest(regresion)
## 
##  studentized Breusch-Pagan test
## 
## data:  regresion
## BP = 0.25949, df = 1, p-value = 0.6105

Teniendo en cuenta que p-valor=0.6105 mayor a 0.05, implica no hay evidencia suficiente de heterocedasticidad.

Normalidad

para comprobar la normalidad vamos a emplear el tes

shapiro.test(resid(regresion))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(regresion)
## W = 0.96613, p-value = 0.01015

Como p-valor=0.06 entonces tiene una tendecia de normalidad, se puede observar esta tendencia en el histograma de residuos.

residuos <- regresion$residuals
hist(residuos, main = "Histograma de los residuos", xlab = "Residuos", col = "lightblue", breaks = 10)

# Grafico de Q-Q

qqnorm(residuos)
qqline(residuos, col = "red")

En la prueba Q-Q plots realizada anteriormente se observa que los puntos Q-Q se acerca a la linea lo que implica un comportamiento normal de los residuos

Conclusiones

En resumen, los resultados obtenidos permiten concluir la veracidad del modelo es aceptable ya que cumple con una correlación mayo a 0,5, además presenta normalidad, independencia de los errores y homocedasticidad. Se recomienda aplicar otro método estadístico que se ajuste mejor al comportamiento de los datos.

Taller 2

Regresión multiples

library(readxl)
datos<-read_excel("C:/Users/Acer/Downloads/taller.xlsx")
colnames(datos)
## [1] "Grosor_inferior" "Grosor_superior" "Altura"          "Peso"           
## [5] "clase_de_cafe"
head(datos)
## # A tibble: 6 × 5
##   Grosor_inferior Grosor_superior Altura  Peso clase_de_cafe
##             <dbl>           <dbl>  <dbl> <dbl> <chr>        
## 1            12.9             5.9   204.   854 bourbon      
## 2            11.5             5.9   177.   670 castillo     
## 3            11.8             5.9   181.   722 castillo     
## 4            11.5             4.5   206.   817 castillo     
## 5            13.8             4.8   240.   653 bourbon      
## 6            13.3             5.9   168.   848 castillo

correlación

cor.test(datos$Grosor_inferior, datos$Peso)
## 
##  Pearson's product-moment correlation
## 
## data:  datos$Grosor_inferior and datos$Peso
## t = 15.015, df = 100, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7610133 0.8837385
## sample estimates:
##       cor 
## 0.8323083
cor.test(datos$Grosor_superior,datos$Peso)
## 
##  Pearson's product-moment correlation
## 
## data:  datos$Grosor_superior and datos$Peso
## t = 0.026542, df = 100, p-value = 0.9789
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1919198  0.1970274
## sample estimates:
##         cor 
## 0.002654172
cor.test(datos$Altura,datos$Peso)
## 
##  Pearson's product-moment correlation
## 
## data:  datos$Altura and datos$Peso
## t = 7.9818, df = 100, p-value = 2.488e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.488630 0.729766
## sample estimates:
##       cor 
## 0.6238251

De acuerdo con la información obtenida se determina lo siguiente: Las variables gruesor inferior y peso tiene una correlacion de 0,83, esta siendo mayor a 0,5 se considera que las variables tiene una correlacion positiva fuerte. Las variables gruesor superior y peso tiene una correlacion apriximada de 0,026,ya que es un valor menor a 0,5 se considera que las variables tiene una correlacion positiva debil. Las variables altura y peso tiene una correlacion apriximada de 0,62, esta siendo mayor a 0,5 se considera que las variables tiene una correlacion positiva fuerte.

modelo de regresión múltiple

modelo<-lm(datos$Peso ~ datos$Grosor_inferior+datos$Grosor_superior+datos$Altura , data=datos)
summary(modelo)
## 
## Call:
## lm(formula = datos$Peso ~ datos$Grosor_inferior + datos$Grosor_superior + 
##     datos$Altura, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -663.28  -83.83  -10.53   55.98  711.70 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -1763.309    157.398 -11.203  < 2e-16 ***
## datos$Grosor_inferior   143.626     11.881  12.089  < 2e-16 ***
## datos$Grosor_superior    31.851     15.796   2.016   0.0465 *  
## datos$Altura              3.928      0.573   6.855 6.39e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 167.9 on 98 degrees of freedom
## Multiple R-squared:  0.8001, Adjusted R-squared:  0.794 
## F-statistic: 130.8 on 3 and 98 DF,  p-value: < 2.2e-16
modelo2<-lm(datos$Peso ~ datos$Altura+datos$Grosor_superior+datos$Grosor_inferior , data=datos)
summary(modelo2)
## 
## Call:
## lm(formula = datos$Peso ~ datos$Altura + datos$Grosor_superior + 
##     datos$Grosor_inferior, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -663.28  -83.83  -10.53   55.98  711.70 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -1763.309    157.398 -11.203  < 2e-16 ***
## datos$Altura              3.928      0.573   6.855 6.39e-10 ***
## datos$Grosor_superior    31.851     15.796   2.016   0.0465 *  
## datos$Grosor_inferior   143.626     11.881  12.089  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 167.9 on 98 degrees of freedom
## Multiple R-squared:  0.8001, Adjusted R-squared:  0.794 
## F-statistic: 130.8 on 3 and 98 DF,  p-value: < 2.2e-16
modelo3<-lm(datos$Peso ~datos$Grosor_inferior+datos$Altura+datos$Grosor_superior , data=datos)
summary(modelo3)
## 
## Call:
## lm(formula = datos$Peso ~ datos$Grosor_inferior + datos$Altura + 
##     datos$Grosor_superior, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -663.28  -83.83  -10.53   55.98  711.70 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -1763.309    157.398 -11.203  < 2e-16 ***
## datos$Grosor_inferior   143.626     11.881  12.089  < 2e-16 ***
## datos$Altura              3.928      0.573   6.855 6.39e-10 ***
## datos$Grosor_superior    31.851     15.796   2.016   0.0465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 167.9 on 98 degrees of freedom
## Multiple R-squared:  0.8001, Adjusted R-squared:  0.794 
## F-statistic: 130.8 on 3 and 98 DF,  p-value: < 2.2e-16

Interpretacion detallada de los coeficientes

coefisientes<-summary(modelo)$coefficients
coefisientes
##                           Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)           -1763.308675 157.3982878 -11.202845 3.036175e-19
## datos$Grosor_inferior   143.626208  11.8809648  12.088766 3.887846e-21
## datos$Grosor_superior    31.851280  15.7958772   2.016430 4.649143e-02
## datos$Altura              3.927913   0.5730046   6.854942 6.393340e-10

Interpretación:

-Intercepto (-2018.169): predice el peso cuando la altura, grosor inferior y grosor superior del tallo son cero. esta parametro no es interpretable en la práctica, pero sirve para definir la recta del modelo.

-grosor inferior (128.934): por cada unidad de medida(cm) del grosor del tallo el peso de este aumenta en 128.9 gramo, si se mantiene constante la altura y el grosor superior del tallo. Este coeficiente es estadisticamente significativo para el ajuste del modelo (p=2.238224e-18).

-Altura (5.475): por cada unidad de medida(cm) de la altura el peso de este aumenta en 5.475 gramo, si se mantiene constante el grosor inferior y el grosor superior del tallo. Este coeficiente es estadisticamente significativo para el ajuste del modelo (p=2.189907e-12).

-Grosor inferior (56.001): por cada unidad de medida(cm) del grosor superior del tallo el peso de este aumenta en 56.001 gramo, si se mantiene constante la altura y el grosor inferior del tallo. Este coeficiente es estadisticamente tambien se considera significativo para el ajuste del modelo (p=8.542353e-04).

Ecuación del modelo ajustados

cat("Peso =", round(coef(modelo)[1], 2), "+",
    round(coef(modelo)[2], 2), "* Gruesor inferior +",
    round(coef(modelo)[3], 2), "* gruesor superior +",
    round(coef(modelo)[4], 5), "* Altura")
## Peso = -1763.31 + 143.63 * Gruesor inferior + 31.85 * gruesor superior + 3.92791 * Altura

Diagnostico del modelo

Grafico de residuos vs valores ajustados

plot(modelo$fitted.values, modelo$residuals,
     main = "Residuos vs Valores Ajustados",
     xlab = "Valores Ajustados", ylab = "Residuos", pch = 19, col = "red")
abline(h = 0, col = "blue", lty = 2)

Analisis la distribución de los punto en la linea cero indica que no se presentan patrones sistemáticos en los residuos, lo cual se considera viable para continuar con el estudio.

Normalidad de los residuos

hist(modelo$residuals, col = "lightblue", main = "Histograma de residuos")

El histograma demuestra una tendensia de los residuos hacia una normalidad.

normal Q-Q plos

qqnorm(modelo$residuals)
qqline(modelo$residuals, col = "red")

interpretación: los residuos siguen una distribución normal en el Q_Q plot ya que se encuentan alineados con la recta de tendecia.

Pruebas complementarias

Homocedasticidad y autocorrelación

library(lmtest)
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 5.6812, df = 3, p-value = 0.1282
dwtest(modelo)
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.8197, p-value = 0.1827
## alternative hypothesis: true autocorrelation is greater than 0

análisis 1.Homocedasticidad:Teniendo en cuenta la prueba de Breusch pragan obtenemos un p-valor=0.6105 mayor a 0.05, implica no hay evidencia suficiente de heterocedasticidad. 2.autorrelación: De acuerdo al estadístico de Durbin-Watson el valor obtenido se aproximó a 2 y el p-valor es mayor a 0.05 lo que indica que no se puede rechazar la hipotesis nula (H0) de independencia de los residuos, por lo tanto, se cumple el supuesto de independencia.

Multicolinealidad

vif_manual <- function(modelo) {
  X <- model.matrix(modelo)[, -1]
  sapply(1:ncol(X), function(i) {
    rsq <- summary(lm(X[, i] ~ X[, -i]))$r.squared
    1 / (1 - rsq)
  })
}

vif_manual(modelo)
## [1] 1.419020 1.526557 1.791990

Interpretación teniendo en cuenta los resultados obtenidos demuestran que VIF se encuentran en un rango de 1 a 2, se considera con una colinealidad moderad; lo cual no implica un problema de multicolinealidad

Anova

anova(modelo)
## Analysis of Variance Table
## 
## Response: datos$Peso
##                       Df  Sum Sq Mean Sq  F value    Pr(>F)    
## datos$Grosor_inferior  1 9580233 9580233 339.6440 < 2.2e-16 ***
## datos$Grosor_superior  1  159609  159609   5.6586   0.01931 *  
## datos$Altura           1 1325439 1325439  46.9902 6.393e-10 ***
## Residuals             98 2764256   28207                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modelo3)
## Analysis of Variance Table
## 
## Response: datos$Peso
##                       Df  Sum Sq Mean Sq F value    Pr(>F)    
## datos$Grosor_inferior  1 9580233 9580233 339.644 < 2.2e-16 ***
## datos$Altura           1 1370360 1370360  48.583 3.704e-10 ***
## datos$Grosor_superior  1  114688  114688   4.066   0.04649 *  
## Residuals             98 2764256   28207                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modelo2)
## Analysis of Variance Table
## 
## Response: datos$Peso
##                       Df  Sum Sq Mean Sq F value    Pr(>F)    
## datos$Altura           1 5381872 5381872 190.801 < 2.2e-16 ***
## datos$Grosor_superior  1 1561331 1561331  55.353 3.893e-11 ***
## datos$Grosor_inferior  1 4122077 4122077 146.138 < 2.2e-16 ***
## Residuals             98 2764256   28207                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova2<- aov(datos$Peso ~ datos$Grosor_inferior*datos$Grosor_superior*datos$Altura , data=datos)
summary(anova2)
##                                                          Df  Sum Sq Mean Sq
## datos$Grosor_inferior                                     1 9580233 9580233
## datos$Grosor_superior                                     1  159609  159609
## datos$Altura                                              1 1325439 1325439
## datos$Grosor_inferior:datos$Grosor_superior               1  116966  116966
## datos$Grosor_inferior:datos$Altura                        1  165887  165887
## datos$Grosor_superior:datos$Altura                        1   80728   80728
## datos$Grosor_inferior:datos$Grosor_superior:datos$Altura  1   12390   12390
## Residuals                                                94 2388285   25407
##                                                          F value   Pr(>F)    
## datos$Grosor_inferior                                    377.066  < 2e-16 ***
## datos$Grosor_superior                                      6.282   0.0139 *  
## datos$Altura                                              52.168 1.32e-10 ***
## datos$Grosor_inferior:datos$Grosor_superior                4.604   0.0345 *  
## datos$Grosor_inferior:datos$Altura                         6.529   0.0122 *  
## datos$Grosor_superior:datos$Altura                         3.177   0.0779 .  
## datos$Grosor_inferior:datos$Grosor_superior:datos$Altura   0.488   0.4867    
## Residuals                                                                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Grosor_inferior la suma de cuadrados y la media cuadrática es 9580233, un estadístico de 339.6440 y un valor-p de p=2.2e-16. esto me indica una relación altamente significativa. el grosor inferior es el más influyente en el modelo Altura la suma de cuadrados es 1325439, un estadístico de 46.9902 y un valor-p de p=6.393e-10. esto me indica una relación altamente significativa. la altura también es muy influyente en el modelo. Grosor_superior la suma de cuadrados es 159609, un estadístico de 5.6586 y un valor-p de p=0.01931. esto me indica una relación significativa. el grosor superior es la variable que menos influye en el modelo residuales la suma de cuadrados es 2764256 y la media cuadrática es 28207. representa la variabilidad no explicada por el modelo.

Conclusión

1 Acuerdo con el modelo obtenido se concluye que las variables altura y grosor inferior del tallo tienen una relación positiva y significativa sobre el peso de estos, es decir, que mientra que aumente estas dos variables el peso del tallo también lo hará. a diferencia de la variables grosor superior del tallo que no muestra tener significacia a la hora de encontrar el peso.

2 En resumen, los resultados obtenidos permiten concluir la veracidad del modelo es aceptable ya que cumple con una correlación mayo a 0,5, además presenta normalidad, independencia de los errores y homocedasticidad.

3 las tres variables son signifivativas en el modelo porque (p<0,05)

Taller 3

Regresión logística.

para nuestro estudio se tendra las siguientes variables. . peso del tallo . clase de café (0=bourbon ; 1=castillo).

estamos interesados en la probalidad de que una tallo de cafe sea de tipo castillo(castillo=1), convertimos dos variables categoricas castillo y bourbon a tipo factor.

#install.packages("tidyverse")
library(tidyverse)
#install.packages("ISLR")
library(ISLR)

datos$clase_de_cafe <- factor(datos$clase_de_cafe, labels = c("bourbon", "castillo"))
datos$clase_de_cafe_num <- ifelse(datos$clase_de_cafe == "castillo", 1, 0)
head(datos)
## # A tibble: 6 × 6
##   Grosor_inferior Grosor_superior Altura  Peso clase_de_cafe clase_de_cafe_num
##             <dbl>           <dbl>  <dbl> <dbl> <fct>                     <dbl>
## 1            12.9             5.9   204.   854 bourbon                       0
## 2            11.5             5.9   177.   670 castillo                      1
## 3            11.8             5.9   181.   722 castillo                      1
## 4            11.5             4.5   206.   817 castillo                      1
## 5            13.8             4.8   240.   653 bourbon                       0
## 6            13.3             5.9   168.   848 castillo                      1

se estudia la relacion que tiene la clase de café respecto al peso del tallo de este. # modelo lineal logistico

modelo_lineal <- lm( datos$clase_de_cafe_num ~ datos$Peso, data = datos)
summary(modelo_lineal)
## 
## Call:
## lm(formula = datos$clase_de_cafe_num ~ datos$Peso, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72303 -0.42967 -0.09859  0.43607  0.78588 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.9330852  0.1358893   6.867 5.64e-10 ***
## datos$Peso  -0.0004774  0.0001264  -3.776  0.00027 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4701 on 100 degrees of freedom
## Multiple R-squared:  0.1248, Adjusted R-squared:  0.1161 
## F-statistic: 14.26 on 1 and 100 DF,  p-value: 0.0002704
library(ggplot2)
ggplot(datos, aes(x = Peso, y = clase_de_cafe_num)) +
  geom_point(aes(color = ""), shape = 1) +
  geom_smooth(method = "lm", color = "gray", se = FALSE) +
  theme_bw() +
  labs(title = "Regresión lineal para clase de café", y = "Probabilidad de Manual") +
  theme(legend.position = "none")
## `geom_smooth()` using formula = 'y ~ x'

Un modelo de regresion lineal no es el adecudo porque los datos de la información no tiende un comportamiento lineal.

predict(modelo_lineal, newdata = data.frame( Peso = 7))
## Warning: 'newdata' had 1 row but variables found have 102 rows
##           1           2           3           4           5           6 
##  0.52538475  0.61322654  0.58840169  0.54304859  0.62134236  0.52824916 
##           7           8           9          10          11          12 
##  0.48480566  0.63757399  0.55593842  0.59794971  0.63566439  0.62802597 
##          13          14          15          16          17          18 
##  0.63327738  0.46618702  0.30625767  0.63327738  0.59842711  0.50963051 
##          19          20          21          22          23          24 
##  0.46618702  0.52538475  0.72302878  0.21077746  0.63900620  0.14728312 
##          25          26          27          28          29          30 
##  0.56214463  0.64043840  0.71109375  0.44040736  0.61083954  0.37595822 
##          31          32          33          34          35          36 
##  0.26902039  0.52061074  0.33728874  0.66526325  0.66287625  0.32344411 
##          37          38          39          40          41          42 
##  0.38741585  0.42513053  0.39696387  0.62038756  0.24467294  0.58887909 
##          43          44          45          46          47          48 
##  0.50819831  0.45281979  0.64664461  0.41605991  0.58219547  0.54161638 
##          49          50          51          52          53          54 
##  0.66192145 -0.08903039  0.60415592  0.37691302  0.07328596  0.37070681 
##          55          56          57          58          59          60 
##  0.26472378  0.41176330  0.57503446  0.41987912  0.25469836  0.45663900 
##          61          62          63          64          65          66 
##  0.32010230  0.43467855  0.18786221  0.42751753  0.43038194  0.50628871 
##          67          68          69          70          71          72 
##  0.36688760  0.54209379  0.64378021  0.41701471  0.47143843  0.56787344 
##          73          74          75          76          77          78 
##  0.57312485 -0.08759819  0.41987912  0.58553728  0.38598364  0.36927461 
##          79          80          81          82          83          84 
##  0.34158535  0.57694406  0.43372375  0.44947798  0.47668984  0.48337346 
##          85          86          87          88          89          90 
##  0.54782260  0.10861364  0.32439891  0.33203733  0.41033110  0.57789886 
##          91          92          93          94          95          96 
##  0.56023503  0.22271249  0.31628309  0.50390170  0.08856279  0.58983389 
##          97          98          99         100         101         102 
##  0.21411927  0.66383105  0.19788763  0.37882263  0.63423219 -0.12483547

Modelo de regresión logistico

modelo_logistico <- glm(datos$clase_de_cafe_num ~ datos$Peso, data = datos, family = binomial)
summary(modelo_logistico)
## 
## Call:
## glm(formula = datos$clase_de_cafe_num ~ datos$Peso, family = binomial, 
##     data = datos)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.1773605  0.7233709   3.010 0.002612 ** 
## datos$Peso  -0.0024224  0.0007291  -3.322 0.000893 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 140.42  on 101  degrees of freedom
## Residual deviance: 126.07  on 100  degrees of freedom
## AIC: 130.07
## 
## Number of Fisher Scoring iterations: 4

La ecuación del modelo logistico es \(logit(P)= 2.177 - 0.0024(Peso)\) donde: .P: probabilidad de que un planta de café sea de castillo. .peso: peso del tallo de café.

El coeficiente estimado para la intersección es el valor esperado de logaritmo de odds de que un tallo de cafe sea bourbon. como es de esperar, los odda son pequeños. si tengo un predictor significativo la devianza tienende a disminuir, lo que sucede en el modelo,entonces considero que hay una contribucion significativa # Grafica

ggplot(datos, aes(x = Peso , y =clase_de_cafe_num )) +
  geom_point(aes( color=as.factor(clase_de_cafe_num)), shape = 1) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  theme_bw() +
  labs(title = "Curva logística con geom_smooth", y = "Probabilidad de Manual") +
  theme(legend.position = "none")
## `geom_smooth()` using formula = 'y ~ x'

##

exp(modelo_logistico$coefficients)
## (Intercept)  datos$Peso 
##   8.8229872   0.9975805

odd de la variable independiente es 0.997, por cada unidad que aumenta la variable peso, el odds de que se presente el evento aumenta 0.997,aproximadamente 1, lo que significa que el aumento es mínimo

predicción y clasificación

Calculamos la probabilidad predicha para cada observación y la clasificamos como “castilla” si la probabilidad es mayor a 0.5.

probabilidades <- predict(modelo_logistico, type = "response")

pred_clase <- ifelse(probabilidades > 0.5, "castillo", "bourbon")

head(data.frame(
  Peso = datos$Peso,
  Probabilidad = round(probabilidades, 3),
  Clasificación = pred_clase
))
##   Peso Probabilidad Clasificación
## 1  854        0.527      castillo
## 2  670        0.635      castillo
## 3  722        0.605      castillo
## 4  817        0.549      castillo
## 5  653        0.645      castillo
## 6  848        0.531      castillo

En la tabla se observa las predicciones maryor de 0.5 de la clase de café castillo

Evaluación de modelo

table(Predicho = pred_clase, Real = datos$clase_de_cafe)
##           Real
## Predicho   bourbon castillo
##   bourbon       36       19
##   castillo      20       27
# Matriz de confusión como tabla
conf <- table(Predicho = pred_clase, Real = datos$clase_de_cafe)

# Convertir a data frame para graficar
conf_df <- as.data.frame(conf)
names(conf_df) <- c("Predicho", "Real", "Frecuencia")

# Gráfico de calor
ggplot(conf_df, aes(x = Real, y = Predicho, fill = Frecuencia)) +
  geom_tile(color = "white") +
  geom_text(aes(label = Frecuencia), size = 6, color = "black") +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(
    title = "Matriz de confusión",
    x = "Clase real",
    y = "Clase predicha"
  ) +
  theme_minimal()

La concordancia de las prediciones indica: .36 plantas no son de cafe de tipo bourbon y efectivamente no lo son. . 20 me indaca que no es de tipo bourbon pero si lo son. . 19 son de de castillo pero en realidad no lo son. . 27 son de castillo y efectivamente lo son.

# Precisión del modelo
mean(pred_clase == datos$clase_de_cafe)
## [1] 0.6176471

La predicion del modelo es de \(61.76%\) considerandoce un aceptable, pero no lo suficiente. ## visualización de el modelo

ggplot(datos, aes(x = Peso, y = as.numeric(clase_de_cafe) - 1)) +
  geom_point(aes(color=as.factor(clase_de_cafe)), size = 3) +
  stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  labs(
    title = "Probabilidad de clase de cafe según peso de los tallos del arbol de café",
    x = "Peso del tallo de café",
    y = "Probabilidad estimada",
    color = "Transmisión"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

conclusiones

.EL modelo muestra que el tallo de la planta de café tiene un efecto negativo y significativo sobre la probabilidad de que clase de cafe sea de tipo castillo.

.A medida que el peso del tallo aumenta la probabilidad de que este sea de tipo castillo disminuye, es decir, que tienen una relación inversamente proporcional.

. De acuerdo con la veracidad arrojada por el modelo se considera aceptable, sin embargo no seria sufuciente para hacer predicciones con rigurocidad, se recomienda hacer un ajuste a la información para el desarrollo del estudio.