Los datos de este fichero se refieren al perímetro, altura y volumen de madera de 31 cerezos negros. Puedes encontrar aquí más información. Las variables son
dbh
: diámetro a la altura del pecho (4 pies y 6 pulgadas, 130 cm) en pies.altura
en pies.volumen
en pies cúbicos.Analizaremos la (posible) relación entre las variables dbh
y volumen
. Se quiere utilizar esta información para estimar el volumen
de madera a partir del dbh
.
Import dataset
que hay en el marco superior derecho de RStudio.read.table()
que se explicó en el vídeo de lectura de datos de la práctica 2. En ese caso, no olvides
cene = read.table(file = "cerezonegro.csv", header = TRUE, sep = ";")
max(cene$dbh)
## [1] 20.6
Comprueba que la lectura ha sido correcta: visualiza la cabecera de la tabla, obtén un resumen de las variables y calcula las dimensiones de la tabla
Para comprobar que la lectura ha sido correcta vamos a ver las primeras lineas de la tabla
head(cene)
## dbh altura volumen
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7
las dimensiones de la tabla y usamos summary
para obtener una exploración inicial de las variables:
dim(cene)
## [1] 31 3
summary(cene)
## dbh altura volumen
## Min. : 8.30 Min. :63 Min. :10.20
## 1st Qu.:11.05 1st Qu.:72 1st Qu.:19.40
## Median :12.90 Median :76 Median :24.20
## Mean :13.25 Mean :76 Mean :30.17
## 3rd Qu.:15.25 3rd Qu.:80 3rd Qu.:37.30
## Max. :20.60 Max. :87 Max. :77.00
Elige las variables explicativa y respuesta. Las variables son diámetro a la altura del pecho y volumen de madera del árbol. Evidentemente, parece más sencillo medir el dimámetro a la altura del pecho, por lo que la elegimos como explicativa. El volumenn será, claro, la variable respuesta.
Visualiza la nube de puntos para determinar de forma visual si tiene sentido utilizar un modelo lineal para relacionar las variables dbh
y volumen
.
Usaremos la función plot
y elegimos para el eje horizontal la variable dbh
, que hemos seleccionado para ser la variable explicativa.
Estudia la existencia de datos atípicos calculando con el boxplot de cada una de las variables.
Dibujamos diagramas de caja (boxplot) de las variables numéricas dbh
y volumen
, con los valores superpuestos. Guardamos cada boxplot en una variable por si hubiera algún valor atípico que quisieramos eliminar
bpdbh = boxplot(cene$dbh)
stripchart(cene$dbh,
method = "jitter", add = TRUE, vertical = TRUE,
pch=19, col= "red", cex=0.5)
bpvolumen = boxplot(cene$volumen)
stripchart(cene$volumen,
method = "jitter", add = TRUE, vertical = TRUE,
pch=19, col= "red", cex=0.5)
Como ves en este último diagrama, hay un valor atípico. El tratamiento de esos valores atípicos es distinto en cada caso y no se pueden establecer reglas generales. A veces un valor atípico se debe simplemente a un error al anotar los datos y, en tal caso, lo más sensato es eliminar esa observación antes de seguir adelante con el análisis. Pero en otras ocasiones un valor atípico apunta a un fenómeno interesante y eliminarlo sería un error.
En este caso, el valor parece plausible y vamos a conservarlo.
Calcula los coeficientes de la recta de regresión. Recuerda que queremos estudiar la posible relación entre dbh
y volumen
. Vamos a construir el modelo de regresión lineal correspondiente usando la función lm
.
(modelo = lm(cene$volumen ~ cene$dbh))
##
## Call:
## lm(formula = cene$volumen ~ cene$dbh)
##
## Coefficients:
## (Intercept) cene$dbh
## -36.943 5.066
fíjate en la sintaxis: a la izquierda está la variable respuesta (dependiente) y a la derecha la explicativa (independiente), como cuando escribes en matemáticas \(y = f(x)\) y extraemos los coeficientes de la recta de regresión con:
modelo$coefficients
## (Intercept) cene$dbh
## -36.943459 5.065856
Como ves, la pendiente de la recta es (aproximadamente) 5.066, mientras que la ordenada en el origen es -36.94.
Recuerda que hemos escrito la ecuación genérica de la recta de regresión así: \[
y = b_0 + b_1 \cdot x
\] siendo \(x\) la variable predictora (en nuestro caso dbh
) e \(y\) la variable respuesta (en nuestro caso volumen
). Como se obtiene:
b0 = modelo$coefficients[1]
b1 = modelo$coefficients[2]
La recta de regresión de volumen
frente a dbh
es, aproximadamente: \[
volumen = -36.94 + 5.066 \cdot dbh
\]
Puedes obtener mucha más información sobre el modelo usando summary
. No te preocupes si en este momento no entiendes mucha de esta información, quedará clara más adelante (en una de las últimas prácticas del curso).
summary(modelo)
##
## Call:
## lm(formula = cene$volumen ~ cene$dbh)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.065 -3.107 0.152 3.495 9.587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.9435 3.3651 -10.98 7.62e-12 ***
## cene$dbh 5.0659 0.2474 20.48 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9331
## F-statistic: 419.4 on 1 and 29 DF, p-value: < 2.2e-16
Vamos a añadir el dibujo de la recta de regresión al diagrama de dispersión usando abline
plot(cene$dbh, cene$volumen, xlab="dbh", ylab="volumen", col="blue", pch=19)
abline(modelo, col="red", lwd=3)
Calcula los valores predichos de la variable volumen
para los valores observados de la variable dbh
.
Para obtener todos los valores predichos por el modelo para los valores observados hay que sustituirlo en la recta. Puedes hacer (no se ejecuta)
modelo$coefficients[1] + modelo$coefficients[2] * cene$dbh
## [1] 5.103149 6.622906 7.636077 16.248033 17.261205 17.767790 18.780962
## [8] 18.780962 19.287547 19.794133 20.300718 20.807304 20.807304 22.327061
## [15] 23.846818 28.406089 28.406089 30.432431 32.458774 32.965360 33.978531
## [22] 34.991702 36.511459 44.110244 45.630001 50.695857 51.709028 53.735371
## [29] 54.241956 54.241956 67.413183
o bien usar la información almacenada en modelo
para hacer la recta de regresión
modelo$fitted.values
## 1 2 3 4 5 6 7 8
## 5.103149 6.622906 7.636077 16.248033 17.261205 17.767790 18.780962 18.780962
## 9 10 11 12 13 14 15 16
## 19.287547 19.794133 20.300718 20.807304 20.807304 22.327061 23.846818 28.406089
## 17 18 19 20 21 22 23 24
## 28.406089 30.432431 32.458774 32.965360 33.978531 34.991702 36.511459 44.110244
## 25 26 27 28 29 30 31
## 45.630001 50.695857 51.709028 53.735371 54.241956 54.241956 67.413183
Calcula el volumen predicho para individuos con diámetro a la altura del pecho de 9, 12.5, 15, 23 pies. Basta con hacer
sort(cene$dbh)
modelo$coefficients[1] + modelo$coefficients[2] * c(9, 12.5, 15, 23)
aunque hay que asegurarse de no extrapolar, es decir, de que los valores del dbh están en el recorrido de la variable.
range(cene$dbh)
## [1] 8.3 20.6
cosa que NO ocurre para el último de ellos. Por tanto, no podemos confiar en el valor de esa predicción.
Calcula e interpreta el coeficiente de correlación. A continuación calculamos el coeficiente de correlación \(r\) de este modelo.
cor(cene$dbh, cene$volumen)
## [1] 0.9671194
Teniendo en cuenta ese valor, el tipo de datos de que se trata y a la vista del gráfico de dispersión, podemos considerar que el ajuste de la recta de regresión a los datos es bueno.
Calcula e interpreta el coeficiente de determinación. Vale
cor(cene$dbh, cene$volumen)^2
## [1] 0.9353199
e indica qué porcentaje de la variabilidad de volumen
se explica a través de dbh
. El resto (de la variabilidad) se explica con otros factores no considerados en el modelo, incluida la variabilidad individual (es decir, el azar). Podrían añadirse variables, es lo que se conoce como regresión lineal múltiple
, una herramienta más avanzada.
{r results='asis', echo=FALSE, comment=NULL} # ExerciseNumber = ExerciseNumber + 1 # cat(paste0(ExerciseString, "**", ExerciseNumber, "**")) #
R
con el nombre iris
.