---
title: "Practica04-2017-18-regresion-no-lineal"
author: "Sección biología"
date: "11 de octubre de 2017"
output: html_document
---

```{r echo=FALSE}
outputType = knitr::opts_knit$get("rmarkdown.pandoc.to")
```


```{r echo=FALSE}
ExerciseString ="### Ejercicio "
ExerciseNumber = 0  
```

## Los datos 
Los datos que aparecen en el fichero [hayas.csv](datos/hayas.csv) se refieren a un estudio realizado sobre desarrollo sostenible de la producción de hayas (Fagus sylvatica) en 86 parcelas de Navarra, parte de los cuales aparecen en la tabla. Sobre cada parcela se han medido las siguientes variables:  

  + `Altura100`: (m / 100 años) Altura dominante a los 100 años.
  + `Edad`: (Años) Edad media de los árboles tipo de la parcela.
  + `Altura dominante`: (metros) la altura media total de los 100   + árboles más gruesos por hectárea.
  + `Altura media`: (metros) altura media que alcanzan los árboles en la parcela.
  + `Diámetro medio`: (centímetros) Diámetro de tronco medio de los árboles de la parcela, medido a 1.3 m por encima del suelo (en la parte más alta de la pendiente).
  + `Densidad`: (pies / Ha) Número de pies de planta (árboles) por unidad de superficie (Hectárea).
  + `Volumen`: (m3 / Ha) Volumen total aprovechable del árbol.
  + `Protección`: SI/NO, según si el espacio está protegido o no.
  + `Calidad`: Alta/Baja, codificada como 1/2.
  
```{r results='asis', echo=FALSE, comment=NULL}
ExerciseNumber = ExerciseNumber + 1
cat(paste0(ExerciseString, ExerciseNumber))
```

**Descarga de los datos.**

+ Selecciona la carpeta de trabajo adecuada para esta sesión con R. 
+ Descarga el fichero y guárdalo en la subcarpeta `datos` de tu carpeta de trabajo**. 
+ Explora el fichero con un editor de texto como el *Bloc de Notas*.
+ Lee ese fichero con R, usando el comando `read.table`  y guárdalo en un data.frame con el nombre `laketrout`. 
+ Comprueba que la lectura ha sido correcta.

```{r, echo = FALSE}
datos = read.table(file = "datos/hayas.csv", sep = ";", header = TRUE)
```

```{r results='asis', echo=FALSE, comment=NULL}
ExerciseNumber = ExerciseNumber + 1
cat(paste0(ExerciseString, ExerciseNumber))
```

Representa la variable `Volumen` en función de la variable  `Diametromedio`. 

```{r, echo = FALSE}
plot(datos$Diametromedio, datos$Volumen)
```
a primera vista:

  * Arriba, a la izquierda, aparece un dato atípico, que procedemos a eliminar
```{r}
(atipico = which(datos$Volumen>600 & datos$Diametromedio<12))
datos = datos[-atipico, ]
```
  
  * La nube de puntos se parece más a una una curva logarítmica que a una recta. Es decir, buscamos `b0` y `b1` tales que 
$$Volumen = b0+b1*\log(Diametromedio)$$

Si defines una nueva variable, llamemosla $z=\log(Diametromedio)$ la expresión anterior es la de una recta
$$Volumen = b0+b1*z$$
Esto sugiere trabajar con la nueva variable `z = log(Diametromedio)` en lugar de con `Diametromedio`. De hecho, si representas `Volumen` en función de `z = log(Diametromedio)`. 

```{r, echo = FALSE}
plot(log(datos$Diametromedio), datos$Volumen)
```
Se observa que los datos aparecen mejor alineados. Puedes comprobar que el coeficiente de correlación es ligeramente mejor cuando usas el logaritmo:
```{r}
modelo = lm(datos$Volumen ~ datos$Diametromedio)
modeloLn = lm(datos$Volumen ~ log(datos$Diametromedio))
signif(sqrt(summary(modelo)$r.squared), digits = 2)
signif(sqrt(summary(modeloLn)$r.squared), digits = 2)
```

De modo que es preferible usar `modeloLn`. 



```{r results='asis', echo=FALSE, comment=NULL}
ExerciseNumber = ExerciseNumber + 1
cat(paste0(ExerciseString, ExerciseNumber))
```

Estima el valor de `volumen` para `Diametromedio`=50. 
Com usaremos la expresión 

$$Volumen = `r signif(modeloLn$coefficients[1], digits= 3)`+`r signif(modelo$coefficients[2], digits= 3)`*\ln(Diametromedio)$$

hay que sustituir en la expresión anterior el logaritmo de `Diametromedio`=50

```{r}
(valor = signif(modeloLn$coefficients[1]+ log(50)*modeloLn$coefficients[2],  digits= 3))

```

Podemos visualizar el resultado

```{r}
par(mfrow = c(1,2))
plot(log(datos$Diametromedio), datos$Volumen)
abline(modeloLn, col = "blue", lwd = 3)
points(log(50), valor, col = "red", pch = 19, cex = 2)

plot(datos$Diametromedio, datos$Volumen)
par(new = TRUE)
D = seq(min(datos$Diametromedio), max(datos$Diametromedio), length.out = 1000)
V =  modeloLn$coefficients[1]+ log(D)*modeloLn$coefficients[2]
plot(D, V, 
     xlim = c(min(datos$Diametromedio), max(datos$Diametromedio)), 
     ylim =c(min(datos$Volumen), max(datos$Volumen)),
     pch = 19, cex = .3, col = "black", xlab="", ylab = "")
points(50, valor, col = "red", pch = 19, cex = 2)
par(mfrow = c(1,1))
```
