Estudiando Agricultura en Taiwán: Curvas con polinomios

Busca la curva con la mínima suma de cuadrados de las desviaciones o la que proporcione el R² más elevado—Experimentación en Agricultura: Curvas de tipo polinómico (pp. 240 – 246).

# Curvas lineal, cuadrática y cúbica 
library(methods)
ggplot(TRONCO, aes(x = ANILLADO, y   = CALLO)) + geom_point()                               + 
  stat_smooth(method = "lm", formula = y~poly(x,1), size = 1, se = FALSE, colour = "black") +  # lineal
  stat_smooth(method = "lm", formula = y~poly(x,2), size = 1, se = FALSE, colour = "blue")  +  # cuadrática
  stat_smooth(method = "lm", formula = y~poly(x,3), size = 1, se = FALSE, colour = "red")   +  # cúbica
  labs(title = "Curvas lineal (negra), cuadrática (azul) y cúbica (roja)" ) +
  ylab("Callo % circunferencia cicatrizada") + xlab("Anchura del anillado (mm)")

En software estadístico 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 POLINÓMICO (pp. 240-246)
#************************************************************

#### Tabla 17.2 (página 243)
input <- ("
           ANILLADO CALLO
           0        100
           3        100
           6        98
           9        89
           12       51
           ")
TRONCO = read.table(textConnection(input), header = TRUE)
TRONCO$ANILLADO = as.numeric(TRONCO$ANILLADO) # de enteros a variables númericas

#### Crear variables cuadrática y cúbica
library(dplyr)

TRONCO =  mutate(TRONCO,
                 ANILLADO2 = ANILLADO*ANILLADO,
                 ANILLADO3 = ANILLADO*ANILLADO*ANILLADO)

library(FSA)

attach(TRONCO); knitr::kable(TRONCO)

ANILLADO	CALLO	ANILLADO2	ANILLADO3
0	100	0	0
3	100	9	27
6	98	36	216
9	89	81	729
12	51	144	1728

# Representación gráfica de datos
# gráfico
plot(ANILLADO, CALLO, pch = 19, 
     xlab = "Anchura del anillado (mm)", 
     ylab = "% circunferencia cicatrizada",
     col  = "blue"); lines(lowess(ANILLADO, CALLO), col = "red")

library(ggplot2)

ggplot(TRONCO, aes(ANILLADO, CALLO)) + geom_point(alpha=2/10, shape=21, fill="blue", colour="black", size=5) + 
  geom_smooth(method="loess", se=TRUE)

# gráfico lineal, cuadrático y cúbico agrupados
library(methods)
ggplot(TRONCO, aes(x = ANILLADO, y   = CALLO)) + geom_point()                               + 
  stat_smooth(method = "lm", formula = y~poly(x,1), size = 1, se = FALSE, colour = "black") +  # lineal
  stat_smooth(method = "lm", formula = y~poly(x,2), size = 1, se = FALSE, colour = "blue")  +  # cuadrática
  stat_smooth(method = "lm", formula = y~poly(x,3), size = 1, se = FALSE, colour = "red")      # cúbica

#### Modelos a comparar
#### Modelo 1: lineal 
modelo.1 <- lm(CALLO ~ ANILLADO,data = TRONCO); summary(modelo.1); anova(modelo.1)

## 
## Call:
## lm(formula = CALLO ~ ANILLADO, data = TRONCO)
## 
## Residuals:
##     1     2     3     4     5 
##  -9.4   1.5  10.4  12.3 -14.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  109.400     10.669  10.254  0.00198 **
## ANILLADO      -3.633      1.452  -2.503  0.08751 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.77 on 3 degrees of freedom
## Multiple R-squared:  0.6761, Adjusted R-squared:  0.5682 
## F-statistic: 6.263 on 1 and 3 DF,  p-value: 0.08751

## Analysis of Variance Table
## 
## Response: CALLO
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## ANILLADO   1 1188.1  1188.1   6.263 0.08751 .
## Residuals  3  569.1   189.7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# gráfico lineal opción 1
plot(ANILLADO, CALLO, pch = 19, 
     xlab = "Anchura del anillado (mm)", 
     ylab = "% circunferencia cicatrizada",
     col  = "blue"); abline(modelo.1, col = "black")

# gráfico lineal opción 2
ggplot(TRONCO, aes(ANILLADO, CALLO)) + geom_point() + 
  geom_smooth(method="lm", formula=y~poly(x,1))

# gráfico lineal opción 3
library(polynom)

library(ggplot2)
VARIABLES <- data.frame("x"=ANILLADO, "y"=CALLO)

my.formula <- y ~ poly(x, 1, raw = TRUE)
p <- ggplot(VARIABLES, aes(x, y)) 
p <- p + geom_point(alpha=2/10, shape=21, fill="blue", colour="black", size=5)
p <- p + geom_smooth(method = "lm", se = FALSE, formula = my.formula, colour = "red")

m <- lm(my.formula, VARIABLES)
my.eq <- as.character(signif(as.polynomial(coef(m)), 3))
label.text <- paste(gsub("x", "~italic(x)", my.eq, fixed = TRUE),
                    paste("italic(R)^2",  format(summary(m)$r.squared, digits = 2), sep = "~`=`~"), sep = "~~~~")

p + annotate(geom = "text", x = 0.0, y = 125, label = label.text, 
             family = "serif", hjust = 0, parse = TRUE, size = 4)

#### Modelo 2: cuadrático
modelo.2 <- lm(CALLO ~ ANILLADO + ANILLADO2, data = TRONCO); summary(modelo.2); anova(modelo.2)

## 
## Call:
## lm(formula = CALLO ~ ANILLADO + ANILLADO2, data = TRONCO)
## 
## Residuals:
##      1      2      3      4      5 
##  2.457 -4.429 -1.457  6.371 -2.943 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  97.5429     5.8406  16.701  0.00357 **
## ANILLADO      4.2714     2.3062   1.852  0.20520   
## ANILLADO2    -0.6587     0.1843  -3.574  0.07014 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.206 on 2 degrees of freedom
## Multiple R-squared:  0.9562, Adjusted R-squared:  0.9123 
## F-statistic: 21.81 on 2 and 2 DF,  p-value: 0.04384

## Analysis of Variance Table
## 
## Response: CALLO
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## ANILLADO   1 1188.10 1188.10  30.848 0.03092 *
## ANILLADO2  1  492.07  492.07  12.776 0.07014 .
## Residuals  2   77.03   38.51                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# gráfico cuadrático opción 1 
plot(ANILLADO, CALLO, pch = 19, 
     xlab = "Anchura del anillado (mm)", 
     ylab = "% circunferencia cicatrizada",
     col  = "blue")

xv <- seq(min(ANILLADO), max(ANILLADO), 0.01)
yv <- predict(modelo.2, list(ANILLADO = xv, ANILLADO2 = xv^2))
lines(xv,yv, col = "blue")

# gráfico cuadrático opción 2 
ggplot(TRONCO, aes(ANILLADO, CALLO)) + geom_point() + 
  geom_smooth(method="lm", formula=y~poly(x,2))

# gráfico cuadrático opción 3 
VARIABLES <- data.frame("x"=ANILLADO, "y"=CALLO)

my.formula <- y ~ poly(x, 2, raw = TRUE)
p <- ggplot(VARIABLES, aes(x, y)) 
p <- p + geom_point(alpha=2/10, shape=21, fill="blue", colour="black", size=5)
p <- p + geom_smooth(method = "lm", se = FALSE, formula = my.formula, colour = "red")

m <- lm(my.formula, VARIABLES)
my.eq <- as.character(signif(as.polynomial(coef(m)), 3))
label.text <- paste(gsub("x", "~italic(x)", my.eq, fixed = TRUE),
                    paste("italic(R)^2",  format(summary(m)$r.squared, digits = 2), sep = "~`=`~"), sep = "~~~~")

p + annotate(geom = "text", x = 0.0, y = 125, label = label.text, 
             family = "serif", hjust = 0, parse = TRUE, size = 4)

#### Modelo 3: cúbico
modelo.3 <- lm(CALLO ~ ANILLADO + ANILLADO2 + ANILLADO3, data = TRONCO); summary(modelo.3); anova(modelo.3)

## 
## Call:
## lm(formula = CALLO ~ ANILLADO + ANILLADO2 + ANILLADO3, data = TRONCO)
## 
## Residuals:
##       1       2       3       4       5 
## -0.2429  0.9714 -1.4571  0.9714 -0.2429 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 100.24286    2.01732  49.691   0.0128 *
## ANILLADO     -2.17857    1.71062  -1.274   0.4238  
## ANILLADO2     0.84127    0.36203   2.324   0.2587  
## ANILLADO3    -0.08333    0.01983  -4.202   0.1487  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.032 on 1 degrees of freedom
## Multiple R-squared:  0.9977, Adjusted R-squared:  0.9906 
## F-statistic: 141.5 on 3 and 1 DF,  p-value: 0.06169

## Analysis of Variance Table
## 
## Response: CALLO
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## ANILLADO   1 1188.10 1188.10 287.775 0.03748 *
## ANILLADO2  1  492.07  492.07 119.187 0.05815 .
## ANILLADO3  1   72.90   72.90  17.657 0.14873  
## Residuals  1    4.13    4.13                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# gráfico cúbico opción 1 
plot(ANILLADO, CALLO, pch = 19, 
     xlab = "Anchura del anillado (mm)", 
     ylab = "% circunferencia cicatrizada",
     col  = "blue")

xv <- seq(min(ANILLADO), max(ANILLADO), 0.01)
yv <- predict(modelo.3, list(ANILLADO = xv, ANILLADO2 = xv^2, ANILLADO3 = xv^3))
lines(xv,yv, col = "red")

# gráfico cúbico opción 2 
ggplot(TRONCO, aes(ANILLADO, CALLO)) + geom_point() + 
  geom_smooth(method="lm", formula=y~poly(x,3))

# gráfico cúbico opción 3 
VARIABLES <- data.frame("x"=ANILLADO, "y"=CALLO)

my.formula <- y ~ poly(x, 3, raw = TRUE)
p <- ggplot(VARIABLES, aes(x, y)) 
p <- p + geom_point(alpha=2/10, shape=21, fill="blue", colour="black", size=5)
p <- p + geom_smooth(method = "lm", se = FALSE, formula = my.formula, colour = "red")

m <- lm(my.formula, VARIABLES)
my.eq <- as.character(signif(as.polynomial(coef(m)), 3))
label.text <- paste(gsub("x", "~italic(x)", my.eq, fixed = TRUE),
                    paste("italic(R)^2",  format(summary(m)$r.squared, digits = 2), sep = "~`=`~"), sep = "~~~~")

p + annotate(geom = "text", x = 0.0, y = 125, label = label.text, 
             family = "serif", hjust = 0, parse = TRUE, size = 4)

#### Comparación de modelos
library(rcompanion)

knitr::kable(compareLM(modelo.1, modelo.2, modelo.3))

Formula

CALLO ~ ANILLADO
CALLO ~ ANILLADO + ANILLADO2
CALLO ~ ANILLADO + ANILLADO2 + ANILLADO3

Rank	Df.res	AIC	AICc	BIC	R.squared	Adj.R.sq	p.value	Shapiro.W	Shapiro.p
2	3	43.86	67.86	42.69	0.6761	0.5682	0.08751	0.9127	0.4839
3	2	35.86	Inf	34.30	0.9562	0.9123	0.04384	0.9339	0.6230
4	1	23.23	-36.77	21.28	0.9977	0.9906	0.06169	0.8810	0.3140

anova(modelo.1, modelo.2, modelo.3)

## Analysis of Variance Table
## 
## Model 1: CALLO ~ ANILLADO
## Model 2: CALLO ~ ANILLADO + ANILLADO2
## Model 3: CALLO ~ ANILLADO + ANILLADO2 + ANILLADO3
##   Res.Df    RSS Df Sum of Sq       F  Pr(>F)  
## 1      3 569.10                               
## 2      2  77.03  1    492.07 119.187 0.05815 .
## 3      1   4.13  1     72.90  17.657 0.14873  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(leaps)

library(HH)

todos <- regsubsets(CALLO ~ ANILLADO + ANILLADO2 + ANILLADO3, data = TRONCO)
knitr::kable(summaryHH(todos))

model	p	rsq	rss	adjr2	cp	bic	stderr
ANILLADO3	2	0.9608872	68.729066	0.9478496	15.647179	-12.98765	4.786407
ANILLADO2-ANILLADO3	3	0.9938397	10.824883	0.9876794	3.621944	-20.61983	2.326465
ANILLADO-ANILLADO2-ANILLADO3	4	0.9976505	4.128571	0.9906019	4.000000	-23.82998	2.031889

En software estadístico SAS

/***********************************************************

EXPERIMENTACIÓN EN AGRICULTURA

Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010

CAPÍTULO 17: REGRESIÓN CURVILÍNEA

CURVAS DE TIPO POLINÓMICO (pp. 240-246)

************************************************************/

DATA TRONCO;

INPUT ANILLADO CALLO;

ANILLADO2 = ANILLADO*ANILLADO;

ANILLADO3 = ANILLADO*ANILLADO*ANILLADO;

LABEL

ANILLADO = "Anchura del anillado (mm) lineal"

ANILLADO2 = "Anchura del anillado (mm) cuadrática"

ANILLADO3 = "Anchura del anillado (mm) cúbica"

CALLO = "% circunferencia cicatrizada";

CARDS;

0 100

3 100

6 98

9 89

12 51

;

ODS HTML;

* Mostrar datos a analizar;

PROC PRINT DATA = TRONCO;

RUN;

*Opción 1;

PROC REG DATA = TRONCO;

MODEL CALLO = ANILLADO;

MODEL CALLO = ANILLADO ANILLADO2;

MODEL CALLO = ANILLADO ANILLADO2 ANILLADO3 / VIF;

RUN;

*Opción 2;

PROC GLM DATA = TRONCO;

MODEL CALLO = ANILLADO ANILLADO*ANILLADO ANILLADO*ANILLADO*ANILLADO / SS1;

RUN;

*Análisis alternativo a STEPWISE;

PROC REG DATA = TRONCO;

MODEL CALLO = ANILLADO ANILLADO2 ANILLADO3 /SELECTION = RSQUARE ADJRSQ CP PRESS MSE SSE;

RUN;

QUIT;

ODS HTML CLOSE;

--------------------------------------------------------------

Archivos adjuntos: capítulo_17b.r, capítulo 17b.sas

Estudiando Agricultura en Taiwán

Páginas

junio 13, 2017

Curvas con polinomios

Formula

No hay comentarios:

Publicar un comentario