Las regresiones curvilíneas (páginas 233 – 239) ajustan mejor los valores de fenómenos naturales: el crecimiento de un patógeno
con el tiempo, por ejemplo.
 |
|
Para comprender estos análisis y sus trasformaciones, recomendamos dar un repaso a los
logaritmos comunes y neperianos, facilitará entender los cambios de base a
exponente.
Códigos
in R
#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010
# CAPÍTULO 17: REGRESIÓN CURVILÍNEA
# CURVAS DE TIPO LOGARÍTMICO Y EXPONENCIAL (pp. 233-239)
#************************************************************
#### Tabla 17.1 (página 235)
MELOCOTON <- data.frame(
HORAS.FRÍO = c(50, 100, 200, 300, 350, 400, 450),
PESO.SECO = c( 3, 4, 6, 8, 14, 30, 71))
attach(MELOCOTON)
#### Transformaciones
LogHORAS.FRÍO = log10(HORAS.FRÍO)
LogPESO.SECO = log10(PESO.SECO)
LnHORAS.FRÍO = log(HORAS.FRÍO)
LnPESO.SECO = log(PESO.SECO)
MELOCOTONERO <- round(data.frame(PESO.SECO, HORAS.FRÍO, LogPESO.SECO, LogHORAS.FRÍO, LnPESO.SECO, LnHORAS.FRÍO),digits = 2)
attach(MELOCOTONERO); knitr::kable(MELOCOTONERO)| PESO.SECO | HORAS.FRÍO | LogPESO.SECO | LogHORAS.FRÍO | LnPESO.SECO | LnHORAS.FRÍO |
|---|---|---|---|---|---|
| 3 | 50 | 0.48 | 1.70 | 1.10 | 3.91 |
| 4 | 100 | 0.60 | 2.00 | 1.39 | 4.61 |
| 6 | 200 | 0.78 | 2.30 | 1.79 | 5.30 |
| 8 | 300 | 0.90 | 2.48 | 2.08 | 5.70 |
| 14 | 350 | 1.15 | 2.54 | 2.64 | 5.86 |
| 30 | 400 | 1.48 | 2.60 | 3.40 | 5.99 |
| 71 | 450 | 1.85 | 2.65 | 4.26 | 6.11 |
#### Representación gráfica de los datos originales
plot(PESO.SECO ~ HORAS.FRÍO, pch = 19, col = "red")
lines(lowess(PESO.SECO ~ HORAS.FRÍO), col = "blue")
#### Regresión lineal con los datos originales
# Fig 17.2. Relación entre el peso seco de yemas de melocotonero y las horas-frío acumuladas durante el invierno
library(ggplot2)
library(ggpmisc)ggplot(MELOCOTONERO, aes(HORAS.FRÍO, PESO.SECO)) +
geom_point(shape = 16, colour = "red", size = 4) +
labs(title = "Fig 17.2. Relación entre el peso seco de yemas y horas-frío")+
ylab("Peso seco de 10 yemas (mg)") + xlab("Horas-frío (horas<7ºC)") +
geom_smooth(method = "loess", size = 1.5, se = T) 
# Regresión sin transformar datos
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = HORAS.FRÍO, y = PESO.SECO)) +
geom_smooth(method = "lm", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# Coeficientes de regresión lineal (página 235)
reg1 <- lm(PESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
reg1$coefficients## (Intercept) HORAS.FRÍO
## -13.1030928 0.1230928summary(reg1)## Call:
## lm(formula = PESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 9.948 4.794 -5.515 -15.825 -15.979 -6.134 28.711
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.10309 14.02402 -0.934 0.3930
## HORAS.FRÍO 0.12309 0.04684 2.628 0.0466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.44 on 5 degrees of freedom
## Multiple R-squared: 0.58, Adjusted R-squared: 0.4961
## F-statistic: 6.906 on 1 and 5 DF, p-value: 0.04665# Gráficos de diagnósticos
par(mfrow=c(1,2))
plot(reg1)
#### Representación gráfica de datos transformados a log10 (página 236)
# Opción 1
reg2 <- lm(LogPESO.SECO ~ LogHORAS.FRÍO, data = MELOCOTONERO)
summary(reg2)## Call:
## lm(formula = LogPESO.SECO ~ LogHORAS.FRÍO, data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.1845 -0.0505 -0.2255 -0.3186 -0.1396 0.1194 0.4302
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7164 0.7749 -2.215 0.0776 .
## LogHORAS.FRÍO 1.1835 0.3301 3.585 0.0158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2856 on 5 degrees of freedom
## Multiple R-squared: 0.7199, Adjusted R-squared: 0.6639
## F-statistic: 12.85 on 1 and 5 DF, p-value: 0.01579# Opción 2
reg2 <- lm(LogPESO.SECO ~ I(LogHORAS.FRÍO), data = MELOCOTONERO) # modelo simple para obtener valores iniciales
summary(reg2)## Call:
## lm(formula = LogPESO.SECO ~ I(LogHORAS.FRÍO), data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.1845 -0.0505 -0.2255 -0.3186 -0.1396 0.1194 0.4302
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7164 0.7749 -2.215 0.0776 .
## I(LogHORAS.FRÍO) 1.1835 0.3301 3.585 0.0158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2856 on 5 degrees of freedom
## Multiple R-squared: 0.7199, Adjusted R-squared: 0.6639
## F-statistic: 12.85 on 1 and 5 DF, p-value: 0.01579# méthod = "lm"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = LogHORAS.FRÍO, y = LogPESO.SECO)) +
geom_smooth(method = "lm", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# método = nls"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = LogHORAS.FRÍO, y = LogPESO.SECO)) +
geom_smooth(method = "nls", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# Gráficos de diagnósticos
plot(reg2)
# Opción 3: Regresión nolineal, empleando "nls"
nolineal <- nls(LogPESO.SECO ~ a + b * LogHORAS.FRÍO, data = MELOCOTONERO, start = list(a = 0, b = 0))
summary(nolineal)##
## Formula: LogPESO.SECO ~ a + b * LogHORAS.FRÍO
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a -1.7164 0.7749 -2.215 0.0776 .
## b 1.1835 0.3301 3.585 0.0158 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2856 on 5 degrees of freedom
##
## Number of iterations to convergence: 1
## Achieved convergence tolerance: 1.895e-08#### Antilogaritmo
# base = b; y = número; x = valor de log
# log_b(x) = y => b^y = x b > 0 y b != 1
# antilog_b(y) = x => b^y = x
#Funcion de antilogaritmo
antilog<-function(lx,base)
{
lbx<-lx/log(exp(1),base=base)
result<-exp(lbx)
result
}
# Antilog: El primer valor representa el intercepto, el segundo la base del logaritmo
antilog(-1.71, 10)## [1] 0.01949845#### Curva de tipo exponencial (página 237)
reg3 <- lm(LogPESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
summary(reg3)##
## Call:
## lm(formula = LogPESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.10072 0.06789 -0.05778 -0.24345 -0.14629 0.03088 0.24804
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2264433 0.1444351 1.568 0.17772
## HORAS.FRÍO 0.0030567 0.0004824 6.336 0.00144 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1796 on 5 degrees of freedom
## Multiple R-squared: 0.8893, Adjusted R-squared: 0.8671
## F-statistic: 40.15 on 1 and 5 DF, p-value: 0.001444# Gráfico método = "lm"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = HORAS.FRÍO, y = LogPESO.SECO)) +
geom_smooth(method = "lm", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# método = "nls"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = HORAS.FRÍO, y = LogPESO.SECO)) +
geom_smooth(method = "nls", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# Gráficos de diagnósticos
plot(reg3)
# Antilog
antilog(0.22, 10) # valor 0.22 y 0.003 obtenidos de los coeficientes de regresión## [1] 1.659587antilog(0.003, 10)## [1] 1.006932#### Curva con logaritmo neperiano en variable dependiente (página 238)
reg4 <- lm(LnPESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
summary(reg4)## Call:
## lm(formula = LnPESO.SECO ~ HORAS.FRÍO, data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 0.22619 0.16474 -0.13814 -0.55103 -0.34247 0.06608 0.57464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.522371 0.332186 1.573 0.17664
## HORAS.FRÍO 0.007029 0.001109 6.335 0.00145 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.413 on 5 degrees of freedom
## Multiple R-squared: 0.8892, Adjusted R-squared: 0.8671
## F-statistic: 40.14 on 1 and 5 DF, p-value: 0.001446# Gráfico método = "lm"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = HORAS.FRÍO, y = LnPESO.SECO)) +
geom_smooth(method = "lm", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# método = "nls"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = HORAS.FRÍO, y = LnPESO.SECO)) +
geom_smooth(method = "nls", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# Gráficos de diagnósticos
plot(reg4)
# Antilog
exp(0.52) # 0.52 es el valor del intercepto## [1] 1.682028#### Curva con logaritmo neperiano en variable independiente (página 239)
reg5 <- lm(PESO.SECO ~ LnHORAS.FRÍO, data = MELOCOTONERO)
summary(reg5)## Call:
## lm(formula = PESO.SECO ~ LnHORAS.FRÍO, data = MELOCOTONERO)
##
## Residuals:
## 1 2 3 4 5 6 7
## 10.674 -1.462 -12.410 -17.916 -14.918 -1.358 37.390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -81.05 56.93 -1.424 0.214
## LnHORAS.FRÍO 18.77 10.53 1.782 0.135
##
## Residual standard error: 21.04 on 5 degrees of freedom
## Multiple R-squared: 0.3885, Adjusted R-squared: 0.2662
## F-statistic: 3.176 on 1 and 5 DF, p-value: 0.1348# Gráfico método = "lm"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = LnHORAS.FRÍO, y = PESO.SECO)) +
geom_smooth(method = "lm", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# método = "nls"
my.formula <- y ~ x
p <- ggplot(data = MELOCOTONERO, aes(x = LnHORAS.FRÍO, y = PESO.SECO)) +
geom_smooth(method = "nls", se=TRUE, color="blue", formula = my.formula) +
stat_poly_eq(formula = my.formula, eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~")),
parse = TRUE) + geom_point() ; p
# Gráficos de diagnósticos
plot(reg5)
# Antilog
exp(-81.0) # -81.0 es el valor intercepto de la ecuación de regresión## [1] 6.639677e-36
----------------------------------------------------
Archivos adjuntos: capítulo 17a.r
Hola, muy buen post, soy agrónomo, no estadístico y estoy aprendiendo a usar R de modo autodidacta, y me surgió el siguiente problema.
ResponderEliminarIngreso las isguientes lineas de código
x<-c(18.26,25.37,16.13,18.96,18.91,19.91)
y<-c(9.11,13.58,9.03,11.48,10.00,9.57)
Reg<-lm(formula=y~x)
anova(Reg)
summary(Reg)
Reg<-lm(formula=y~x-1)
anova(Reg)
summary(Reg)
No obstante, el R^2 y el p-value de la salidas del modelo sin intercepto no convence, aún más cuando lo comparo con la salida de la regresión con intercepto.
Si pudera ayudarme, yo le agradecería mucho.
Saludos!
Hola Víctor,
Eliminar¡Muy buena observación! El R^2 no se calculó correctamente debido a la ausencia de un intercepto. Los resultados que obtuviste con el software son correctos. La diferencia se produce principalmente en el Total de Suma de Cuadrados (SST) de la regresión sin intercepto. Con intercepto el SST es 15.64588; no obstante, sin intercepto el número sube a 672.3247. Te envío la muestra de los cálculos.
#### datos
x<-c(18.26,25.37,16.13,18.96,18.91,19.91)
y<-c(9.11,13.58,9.03,11.48,10.00,9.57)
#### regresiones con (int.reg) y sin intercepto (noint.reg)
int.reg <- lm(y~x) ; anova(int.reg) ; summary(int.reg)
noint.reg <- lm(y~x+0); anova(noint.reg); summary(noint.reg)
#### para comparación
# incluye intercepto
y_ajustados <- int.reg$fitted.values
SSE <- sum((y_ajustados - y)^2); SSE
SST <- sum((y - mean(y))^2); SST
1-SSE/SST # R^2 con intercepto
# sin intercepto
y_ajustados2 <- noint.reg$fitted.values
SSE2 <- sum((y_ajustados2 - y)^2); SSE2 # SSE2 es solo ligeramente mas alto que SSE..
SST2 <- sum((y - 0)^2) ; SST2 # !!! la diferencia esta aqui !!!
1-SSE2/SST2 # R^2 sin intercepto
Encontraras explicaciones teóricas en estos sitios:
https://econintec.files.wordpress.com/2008/06/regresion-a-traves-del-origen.pdf
https://fvela.files.wordpress.com/2010/06/apunte_extensiones-de-regresion.pdf
Y, ejemplos de codigos en los dos siguientes:
https://stats.stackexchange.com/questions/26176/removal-of-statistically-significant-intercept-term-increases-r2-in-linear-mo
https://stackoverflow.com/questions/20333600/why-does-summary-overestimate-the-r-squared-with-a-no-intercept-model-formula
Por su tiempo y explicación ¡Muchas Gracias! Ahora seré un seguidor del blog!!! muy interesante y delicioso el tema de las pithayas!
EliminarSolo otra duda, que no me queda claro con al leer las páginas que me recomendaste.
Desde la perspectiva de la interpretación de resultados.
¿Tiene sentido interpretar el R^2 y p-value, a pesar de estar sobreestimados en el modelo sin intercepto?
¡Gracias!
¡Muy amable!
Tu pregunta tiene implícita solo una respuesta, lo cual no es cierto. Es como si me preguntaras ¿Tiene sentido medir la temperatura únicamente empleando el punto de congelación y ebullición del agua, en grados Celsius? Eso no es exacto, pues también existen otras maneras de cuantificar los cambios de temperatura, además, ya conoces que hay temperaturas más bajas y más altas que el punto de congelación y ebullición; igual, hay otras unidades de medida: Fahrenheit y Kelvin.
EliminarAsí que si desconoces a priori que E [Y] = 0 cuando x = 0 (o no tienes seguridad si la regresión atraviesa el cero en la gráfica x,y = 0,0), entonces probablemente no deberías estar ajustando una línea a través del origen.
Cuando la línea de regresión pasa por el origen, es decir, β0 = 0, entonces R2 no se puede usar para juzgar el modelo ajustado. Para comparar modelos de regresión, te recomiendo analizarlos con otros criterios, como el Criterio de Información de Akaike (AIC) y el Mallows's Cp — ambos indican que variables debes emplear o descartar en tus modelos, y eso incluye el intercepto. Elije el modelo que ofrezca los resultados con menores valores.
En tu ejemplo es la regresión sin intercepto. Aunque la diferencia es muy mínima, con intercepto el AIC=19.812103 y sin intercepto AIC=17.893030. Yo no me complicaría, elegiría cualquiera de los dos. Sin embargo, optaría por la regresión CON INTERCEPTO si estoy seguro que ES IMPOSIBLE que existencia algo que tengo valor cero en la vida real. En el siguiente paper, encontraras un ejemplo, en donde la altura (y) depende del ancho del huevo (x); en él, el autor dice que es imposible que exista un huevo con nada de ancho (ver paginas 78-79): https://online.stat.psu.edu/~ajw13/stat501/SpecialTopics/Reg_thru_origin.pdf
Cuéntame, ¿Qué libro de texto lees? Pregunto porque posiblemente necesites uno que trate modelos lineales de manera más avanzada; fórmulas que emplean el alfabeto griego, con notaciones matemáticas, sumatorias, sigmas,... son ineludibles. Saludos, Edgardo
#### datos
x<-c(18.26,25.37,16.13,18.96,18.91,19.91)
y<-c(9.11,13.58,9.03,11.48,10.00,9.57)
#### regresión con intercepto (int.reg)
# Computo de AIC y Cp de Mallows para regresion con intercepto
int.reg <- lm(y~x) ; anova(int.reg) ; summary(int.reg);
print(sprintf('con intercept AIC=%f',AIC(int.reg)))
MSE=(summary(int.reg)$sigma)^2 # guardar MSE para el modelo completo
intercepto <- lm(y~1) # regresión incluye solo la constante del intercepto
Cp.int.reg <- step(intercepto, scope = list(upper = int.reg), scale = MSE); Cp.int.reg
#### regresión sin intercepto (noint.reg)
# Computo de AIC y Cp de Mallows para regresion sin intercepto
noint.reg <- lm(y~x+0); anova(noint.reg); summary(noint.reg);
print(sprintf('sin intercept AIC=%f',AIC(noint.reg)))
MSE=(summary(noint.reg)$sigma)^2
intercepto <- lm(y~-1) # regresión no incluye intercepto
Cp.noint.reg <- step(intercepto, scope = list(upper = noint.reg), scale = MSE); Cp.noint.reg
Este comentario ha sido eliminado por el autor.
EliminarEste comentario ha sido eliminado por el autor.
EliminarGracias Edgardo, muy amable por tu respuesta.
ResponderEliminarEn cuanto a libro: estoy leyendo The R Book de Crawley y llegué a ese libro al leer la introducción a la estadística con R que está en el CRAN castellano.
Felices fiestas de fin de año.
¡Saludos!
Excelente texto. Tienes muy buena referencia. Igualmente, disfruta las festividades de esta época. E.
Eliminar