#Introducción

En este práctico analizaremos la relación entre la profundidad del suelo y el contenido de nitrógeno (%), usando herramientas de correlación y regresión lineal simple en R.

Carga de Paquetes necesarios

*todas las librerias ya están instaladas, si que solo las cargamos.

library(tidyverse)       # Manipulación de datos y gráficos
library(readxl)          # Lectura de archivos Excel
library(PerformanceAnalytics) # Gráficos de correlación
library(lmtest)          # Test de supuestos en regresión

A. Exploración de Datos

Cargar la base de datos

NITROGENO <- read_excel("NITROGENO.xlsx")

#Explorar la estructura de los datos

glimpse(NITROGENO)
## Rows: 12
## Columns: 3
## $ REPETICION  <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4
## $ PROFUNDIDAD <dbl> 20, 20, 20, 20, 40, 40, 40, 40, 60, 60, 60, 60
## $ NITROGENO   <dbl> 0.09, 0.08, 0.13, 0.11, 0.06, 0.08, 0.09, 0.09, 0.04, 0.05…
summary(NITROGENO)
##    REPETICION    PROFUNDIDAD   NITROGENO      
##  Min.   :1.00   Min.   :20   Min.   :0.04000  
##  1st Qu.:1.75   1st Qu.:20   1st Qu.:0.05750  
##  Median :2.50   Median :40   Median :0.08000  
##  Mean   :2.50   Mean   :40   Mean   :0.07667  
##  3rd Qu.:3.25   3rd Qu.:60   3rd Qu.:0.09000  
##  Max.   :4.00   Max.   :60   Max.   :0.13000

Los datos son correctos, según lo que pude ver de la base de datos “Nitrogeno.xlsx”

Graficar la relación entre profundidad y nitrógeno %

ggplot(NITROGENO, aes(x = PROFUNDIDAD, y = NITROGENO)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "blue") +
  labs(x = "Profundidad (cm)", y = "Nitrógeno (%)") +
  theme_minimal()

Los puntos marcan las muestras que tenemos, agregué una linea de regresión lineal (o de mínimos cuadrados) porque creo que es más sencillo interpretar el resto del practyico de esta forma.

B. Análisis de Correlación

Calcular coeficiente de correlación de Pearson

cor(NITROGENO$PROFUNDIDAD, NITROGENO$NITROGENO, method = "pearson")
## [1] -0.8453206

El coeficiente es alto y de signo negativo, lo que demuestra que hay una correlación fuerte entre las variables. Sin embargo, la correlación es negativa, a mayor profundidad, menor contenido de nitrógeno porcentual. Está de acuerdo a lo que vimos en el gráfico anterior.-

Prueba de hipótesis para ρ

cor.test(NITROGENO$PROFUNDIDAD, NITROGENO$NITROGENO, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  NITROGENO$PROFUNDIDAD and NITROGENO$NITROGENO
## t = -5.0034, df = 10, p-value = 0.0005346
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9556210 -0.5271611
## sample estimates:
##        cor 
## -0.8453206

El valor de p para esta prueba (0.0005346) demuestra que es muy poco probable que estos resultados sean solo causados por la toma azarosa de muestras, lo que sugiere que podemos descartar la hipótesis nula (que en este caso es que no hay correlación)

C. Ajuste del Modelo de Regresión Lineal Simple Ajustar modelo

MODELO <- lm(NITROGENO ~ PROFUNDIDAD, data = NITROGENO)
summary(MODELO)
## 
## Call:
## lm(formula = NITROGENO ~ PROFUNDIDAD, data = NITROGENO)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.024167 -0.010417  0.002083  0.011458  0.025833 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1316667  0.0118732  11.089 6.11e-07 ***
## PROFUNDIDAD -0.0013750  0.0002748  -5.003 0.000535 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01555 on 10 degrees of freedom
## Multiple R-squared:  0.7146, Adjusted R-squared:  0.686 
## F-statistic: 25.03 on 1 and 10 DF,  p-value: 0.0005346

(Intercept) indica el valor de nitrógeno porcentual cuando la profundidad es de 0 cm En este caso 0.13%. No es un valor de real, pero sirve de referencia.

Pendiente (–0.001375) indica que por cada cm de profundidad que ganamos, perdemos 0.001375% en contenido porcentual de Nitrógeno.

La bondad de ajuste o coeficiente de Pearson de R² = 0.715 (71.5%) indica que porcentaje de un valor de contenido de nitrógeno puede ser predicho por el modelo, si se cuenta con la profundidad en cm.

Gráfico con la recta ajustada

ggplot(NITROGENO, aes(x = PROFUNDIDAD, y = NITROGENO)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "blue") +
  labs(x = "Profundidad (cm)", y = "Nitrógeno (%)") +
  theme_minimal()

D. Diagnóstico del Modelo Gráficos de residuos

par(mfrow=c(2,2))
plot(MODELO)

par(mfrow=c(1,1))

no parecen haber discrepancias o comportanmientos extraños en ninguno de los gráficos realizados. esto quiere dcir que podemos confiar en lo que nos diga el modelo.

Normalidad de residuos

shapiro.test(residuals(MODELO))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(MODELO)
## W = 0.97087, p-value = 0.9197

Los residuos son normales, el p mayor que 0.9 quiere decir que es muy probablñe que los residuos tengan una distribución normal. Podemos trabajar con los datos tal y como son, sin transformarlos.

Homocedasticidad (test de Breusch–Pagan)

bptest(MODELO)
## 
##  studentized Breusch-Pagan test
## 
## data:  MODELO
## BP = 4.2387, df = 1, p-value = 0.03951

Hay heterocedasticidad en los datos, hay una profundidd, al menos, en que la variabilidad de los datos no es la misma que en las demas. En este caso, la profuindidad media tiene menor variabilidad. Los datos en conjunto dan un modelo aceptable así que, aunque podría usar una transformación para mejorar mi ajuste, no aprece necesario.

# Modelo con transformación logarítmica
MODELO_LOG <- lm(log(NITROGENO) ~ PROFUNDIDAD, data = NITROGENO)

# Gráfico con recta ajustada en escala original
ggplot(NITROGENO, aes(x = PROFUNDIDAD, y = log(NITROGENO))) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  labs(x = "Profundidad (cm)", 
       y = "log(Nitrógeno %)",
       title = "Relación transformada: log(N) vs Profundidad") +
  theme_minimal()

# Test de normalidad de residuos
shapiro.test(residuals(MODELO_LOG))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(MODELO_LOG)
## W = 0.86346, p-value = 0.05402
# Test de homocedasticidad (Breusch–Pagan)
bptest(MODELO_LOG)
## 
##  studentized Breusch-Pagan test
## 
## data:  MODELO_LOG
## BP = 0.15063, df = 1, p-value = 0.6979

Apliqué una transformación por logaritmos,pero sin mejoras evidentes. Probablemente porque es la profundidad intermedia la que presenta menor variabilidad.

Predicción Predicción para una profundidad de 50 cm

nuevo <- data.frame(PROFUNDIDAD = 50)
predict(MODELO, newdata = nuevo, interval = "confidence")
##          fit        lwr        upr
## 1 0.06291667 0.05119171 0.07464163

el contenido deberia ser de aproximadamente 0.063%, con un mínimo de 0.051% y un máximo de 0.075%