Parcial.3.3
Luis Fernando Noguera
5/28/2020
##Cargar los Datos.
## Primero vamos a cargar la base de datos de todos los puntos para posterior emprezar a resolver punto a punto.
load("/Users/luisnoguera/Desktop/Parcial.3/gala.RData")
datos_plantas = gala
########################################################
load("/Users/luisnoguera/Desktop/Parcial.3/cafe.Rdata")
datos_cafe = cafe_calidad
########################################################
load("/Users/luisnoguera/Desktop/Parcial.3/moluscos.RData")
datos_moluscos = BD_moluscosPUNTO 1: Regresión Lineal
Primero vamos a cargar la base de datos conocida como GALA, la cual nos brinda datos del número de especies de plantas en 30 islas de Galápagos (Faraway, 2005). En donde podemos encontrar el número de especies, clasificadas con las variables Endemic (número de especies endémicas), y estudiadas con variables como el Area (km2), Elevation (altura máxima, m), Nearest (distancia a la isla más cercana, km), Scruz (distancia hasta la isla Santa Cruz, km), Adjacent (área de la isla adyacente, km2). Suponiendo que las 30 observaciones se comportan como una muestra aleatoria, vamos a estimar un modelo de regresión múltiple para estudiar qué variables nos permiten explicar mejor el número de especies.
Realizaremos un analisis bivariado el cual nos brinda informacion necesaria para ver la fuerza con la que estan relaciodas dos variables y la direccion de dicha fuerza, en pocas palabras me permite ver que tanto o poco estan relacionadas unas con otras y si la fuerza es directa o inversa, para asi plantear un modelo de regresion lineal con variables que puedan dar estimaciones acertadas.
Para realizar lo anteriormente descrito realizaremos como primera instacia un acercamiento gráfico a los datos mediante un diagrama de dispersión, para observar las características importantes y la relación entre los mismos, posterior revisaremos el coeficiente de correlacion de las variables y de esta manera podremos determinar el modelo final de regresion multiple.
plot(datos_plantas$Endemics,datos_plantas$Species)
plot(datos_plantas$Area,datos_plantas$Species)
plot(datos_plantas$Elevation,datos_plantas$Species)
plot(datos_plantas$Nearest,datos_plantas$Species)
plot(datos_plantas$Scruz,datos_plantas$Species)
plot(datos_plantas$Adjacent,datos_plantas$Species)
## especies es la variable que quiero explicar Con los resultados graficos podemos ver que el numero de especies se explica bastante bien en terminos de Endemics y elevation. Sin embargo vamos a realizar el analisis de correlacion entre variables para obtener resultados mas concretos.
cor(datos_plantas,datos_plantas$Species,use = "complete.obs")## [,1]
## Species 1.00000000
## Endemics 0.97087652
## Area 0.61784307
## Elevation 0.73848666
## Nearest -0.01409407
## Scruz -0.17114244
## Adjacent 0.02616635Como habiamos visto anteriormente en los graficos el numero de especies se explican bastante bien en terminos de Endemics, Area y Elevation. Por otro lado vemos como el resto de variables son menos importantes para el planteamiento de mi modelo de regresion multiple.
A continuacion realizaremos el modelo de regresion multiple con las variables anteriormente mensionadas, de las que obtuvieron los siguientes resultados.
mod1 = lm(datos_plantas$Species~datos_plantas$Endemics+datos_plantas$Area+datos_plantas$Elevation)
summary(mod1)##
## Call:
## lm(formula = datos_plantas$Species ~ datos_plantas$Endemics +
## datos_plantas$Area + datos_plantas$Elevation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.691 -10.530 2.387 10.529 72.723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.891237 7.569210 -2.099 0.0456 *
## datos_plantas$Endemics 4.331791 0.304686 14.217 8.97e-14 ***
## datos_plantas$Area 0.012669 0.008936 1.418 0.1681
## datos_plantas$Elevation -0.041439 0.023653 -1.752 0.0916 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.29 on 26 degrees of freedom
## Multiple R-squared: 0.9492, Adjusted R-squared: 0.9433
## F-statistic: 161.8 on 3 and 26 DF, p-value: < 2.2e-16Después de realizar la estimación de la Species en función de Endemics, Area y Elevation se obtuvieron los siguientes resultados: p-value: < 2.2e-16 R2 : 0,95 El valor p es menor a 0.05 y muy cercano a 0, lo que significa que rechazamos hipotesis nula y aceptamos hipotesis alterna, en otras palablas intentar hacer la correlación tiene sentido, es decir existe una relación entre las variables Species, Endemics, Area y Elevation.
El R2 permite concluir que el modelo explica o tiene una precisión de un 95% a la hora de simular la Species de la planta en función de Endemics, Area y Elevation, o visto desde otro lado, el modelo tiene un error del 5%.
Finalmente realizaremos una estimacion con el modelo anteriormente calculado.
datos_reales = t(datos_plantas$Species)
datos_est = data.frame(predict(mod1,datos_plantas))
colnames(datos_est)[1] = "Numero de especies estimado"
errores = c(1:30)
datos_est## Numero de especies estimado
## Baltra 69.719986
## Bartolome 70.575249
## Caldwell -7.617234
## Champion 21.189962
## Coamano -14.749605
## Daphne.Major 26.831547
## Daphne.Minor -19.744038
## Darwin 7.499091
## Eden -1.505852
## Enderby -11.866527
## Espanola 89.268685
## Fernandina 81.850429
## Gardner1 55.725928
## Gardner2 -7.960813
## Genovesa 63.483255
## Isabela 358.059788
## Marchena 71.166995
## Onslow -8.263500
## Pinta 112.941626
## Pinzon 108.306284
## Las.Plazas 19.202543
## Rabida 98.916385
## SanCristobal 242.993682
## SanSalvador 304.691338
## SantaCruz 371.276645
## SantaFe 94.971327
## SantaMaria 275.974100
## Seymour 47.349219
## Tortuga 11.071173
## Wolf 25.642331datos_est = t(datos_est)
for (i in c(1:length(datos_reales))) {
errores[i] = datos_reales[i] - datos_est[i]
}
plot(errores, col = 'blue')
plot(datos_reales, datos_est, col = 'blue')
Note que cuando realizamos la prediccion con nuestro modelo obtenemos muy buenos resultados, con errores relativamente pequeños, por lo tanto el modelo es bastante bueno, sin embargo, en el calculo de numero de especies de algunas islas, el error es bastante grande. la grafica tambien muestra como se comportan los errores, y el claro que algunos estan muy lejos.
PUNTO 2: Modelo de Diseño de Experimentos
Para realizar la practica de modelo de diseño de experiemntos vamos a cargar la base de datos moluscos la cual contiene dos tipos de moluscos A y B, los cuales fueron sometidos a tres concentraciones distintas de agua de mar (100%, 75% y 50%), con lo que se llevo acabo la el consumo de oxígeno midiendo la proporción de O2 por unidad de peso seco del molusco.
Primero se realizara un análisis exploratorio que nos permita conocer como es el consumo de oxígeno en las distintas concentraciones de agua de mar y ademas si estas conclusiones son las mismas para cada tipo de molusco.
require(ggplot2)
require(plotly)
consumo = datos_moluscos$cons_o
datos_moluscos$c_agua=as.character(datos_moluscos$c_agua)
cagua=datos_moluscos$c_agua
Grafico1 = ggplot(datos_moluscos,aes(x=consumo,fill=cagua))+geom_histogram()+facet_grid(~datos_moluscos$molusco)+theme_bw()
ggplotly(Grafico1)revisando los graficos anteriores, nos damos cuenta que el mayor cosumo de oxigeno para los moluscos tipo A es a una concentracion de agua de mar del 50%, ademas que para una concentracion del 100% nos encontramos que su consumo se encuentra entre 5 y 10 unidades de O2. Por otro lado el mayor cosumo de oxigeno para los moluscos tipo para los moluscos tipo B es a una concentracion de agua de mar del 50%, ademas que para una concentacion de agua de mar del 75% todos consumen menos de 10 unidades de O2.
boxplot(consumo~datos_moluscos$molusco,col = "blue")
Note que el consumo de los moluscos tipo A es apenas mayor a los consumos de los moluscos tipo B, sin embargo podemos decir que el consumo de oxigeno en cualquiera de las concentraciones, es menor para los moluscos tipo B.
Ahora realizaremos la estimacion de un modelo de diseño de experimentos el cual permita evaluar el efecto de la concentración de agua de mar y los tipos de molusco sobre el consumo de oxigeno.
mod2=lm(datos_moluscos$cons_o~datos_moluscos$c_agua+datos_moluscos$molusco)
summary(mod2)##
## Call:
## lm(formula = datos_moluscos$cons_o ~ datos_moluscos$c_agua +
## datos_moluscos$molusco)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1750 -1.9877 -0.7019 2.1244 6.1450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.3669 0.8521 10.993 3.33e-14 ***
## datos_moluscos$c_agua50 3.5794 1.0436 3.430 0.00132 **
## datos_moluscos$c_agua75 -1.6788 1.0436 -1.609 0.11486
## datos_moluscos$moluscoB -1.3913 0.8521 -1.633 0.10966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.952 on 44 degrees of freedom
## Multiple R-squared: 0.3986, Adjusted R-squared: 0.3575
## F-statistic: 9.719 on 3 and 44 DF, p-value: 4.866e-05Ahora interpretaremos los coeficientes del modelo, donde podemos apreciar que el valor p es 4.866e-05, osea lo suficientemente pequeño como para aaceptar el modelo, sin embargo el modelo no explica muy bien el consumo de oxigeno de los moluscos.
PUNTO 3: Analisis de Componentes Principales
Para el analisis de componentes principales se trabajara con una base de datos llamada cafe la cual contiene la seleccion un total de 10 cafés que fueron evaluados en cuanto a características de calidad por un grupo de catadores expertos quienes evaluaron en una escala de cero a diez, algunos factores como la intensidad y el aroma del café. Estos datos se encuentran en la base “cafe_calidad” que se le comparte dentro del archivo “cafe.Rdata”. Con base en estos datos realice un análisis de componentes principales con el cual logre identificar las correlaciones entre las variables (factores evaluados) y similitudes entre tipos de café.
Antes de realizar la descomposición espectral de la matriz de varianzas y covarianzas, veamos la matriz de correlacion que existe entre los datos.
## Intensidad Aroma Cuerpo Acidez Amargo Astringencia
## Intensidad 1.00 0.83 0.84 0.87 0.70 0.78
## Aroma 0.83 1.00 0.86 0.72 0.71 0.66
## Cuerpo 0.84 0.86 1.00 0.67 0.66 0.62
## Acidez 0.87 0.72 0.67 1.00 0.67 0.61
## Amargo 0.70 0.71 0.66 0.67 1.00 0.56
## Astringencia 0.78 0.66 0.62 0.61 0.56 1.00Posteriormente, veamos la descomposición espectral de la matriz de varianzas y covarianzas, para realizar la interpretacion correspondiente de los coeficientes.
## Duality diagramm
## class: pca dudi
## $call: dudi.pca(df = datos_cafe, scannf = FALSE, nf = 3)
##
## $nf: 3 axis-components saved
## $rank: 6
## eigen values: 4.601 0.4694 0.3845 0.3446 0.1449 ...
## vector length mode content
## 1 $cw 6 numeric column weights
## 2 $lw 10 numeric row weights
## 3 $eig 6 numeric eigen values
##
## data.frame nrow ncol content
## 1 $tab 10 6 modified array
## 2 $li 10 3 row coordinates
## 3 $l1 10 3 row normed scores
## 4 $co 6 3 column coordinates
## 5 $c1 6 3 column normed scores
## other elements: cent norm## Inertia information:
## Call: inertia.dudi(x = acp_cafe)
##
## Decomposition of total inertia:
## inertia cum cum(%)
## Ax1 4.60148 4.601 76.69
## Ax2 0.46937 5.071 84.51
## Ax3 0.38451 5.455 90.92
## Ax4 0.34462 5.800 96.67
## Ax5 0.14488 5.945 99.08
## Ax6 0.05514 6.000 100.00## inertia cum cum(%)
## Ax1 4.60147661 4.601477 76.69128
## Ax2 0.46937103 5.070848 84.51413
## Ax3 0.38451212 5.455360 90.92266
## Ax4 0.34461612 5.799976 96.66626
## Ax5 0.14487978 5.944856 99.08093
## Ax6 0.05514434 6.000000 100.00000Escojo solo las tres primeras componentes principales, pues con las tres primeras, el acumulado es igual a 90.9 lo que significa que los tres primeros valores principales representa en un 91% el total de los datos.
El circulo de correlaciones, nos esta diciendo que la acidez y aroma estan fuertemente correlacionadas y ademas en un sentido directo, note que ademas todas estan en el cuadrante uno y cuatro del circulo, por lo que se encuentran relacionadas de manera directa, ninguna esta opuesta por lo tanto ninguna se relaciona fuerte de manera inversa.
Note que, podemos ver la buena calificacion que tienen los excelsos tanto claro como oscuros,por su aroma, acidez, astringencia, intensidad y amargo, sin embargo podriamos decir que el excelso claro tiene mas astringencia y que el excelso oscuro tiene mas amargo. Por otro lado veamos la relacion de cada tipo de cafe.
tipo=c("Cl","M","C","M","C","O","M","C","M","Cl")
tipo=as.factor(tipo)
s.class(acp_cafe$li,tipo)
Ahora veamos la base de datos completa, en donde vamos a tener en cuenta tanto las clasificaciones quimicas como las de calidad, para realizar la correspondiente descomposición espectral de la matriz de varianzas y covarianzas.
cafe_quimica[,7:12] =datos_cafe
cafe = cafe_quimica
cafe_quimica=cafe[,1:6]
###### Matriz de Correlaciones ######
data.frame(round(cor(cafe),2))## Color DA pH AcidezT Cafeina AcidosCl Intensidad Aroma Cuerpo
## Color 1.00 0.79 -0.71 0.11 -0.32 0.01 -0.46 -0.76 -0.68
## DA 0.79 1.00 -0.27 -0.34 -0.63 -0.41 -0.74 -0.81 -0.78
## pH -0.71 -0.27 1.00 -0.67 -0.29 -0.60 -0.22 0.17 0.03
## AcidezT 0.11 -0.34 -0.67 1.00 0.88 0.98 0.78 0.47 0.48
## Cafeina -0.32 -0.63 -0.29 0.88 1.00 0.94 0.88 0.79 0.67
## AcidosCl 0.01 -0.41 -0.60 0.98 0.94 1.00 0.81 0.57 0.53
## Intensidad -0.46 -0.74 -0.22 0.78 0.88 0.81 1.00 0.83 0.84
## Aroma -0.76 -0.81 0.17 0.47 0.79 0.57 0.83 1.00 0.86
## Cuerpo -0.68 -0.78 0.03 0.48 0.67 0.53 0.84 0.86 1.00
## Acidez -0.25 -0.61 -0.45 0.82 0.89 0.88 0.87 0.72 0.67
## Amargo -0.63 -0.80 0.07 0.33 0.53 0.38 0.70 0.71 0.66
## Astringencia -0.51 -0.63 0.07 0.60 0.76 0.63 0.78 0.66 0.62
## Acidez Amargo Astringencia
## Color -0.25 -0.63 -0.51
## DA -0.61 -0.80 -0.63
## pH -0.45 0.07 0.07
## AcidezT 0.82 0.33 0.60
## Cafeina 0.89 0.53 0.76
## AcidosCl 0.88 0.38 0.63
## Intensidad 0.87 0.70 0.78
## Aroma 0.72 0.71 0.66
## Cuerpo 0.67 0.66 0.62
## Acidez 1.00 0.67 0.61
## Amargo 0.67 1.00 0.56
## Astringencia 0.61 0.56 1.00acp_cafe=dudi.pca(cafe,scannf= FALSE, nf=5)
acp_inertia=inertia.dudi(acp_cafe)
acp_inertia$tot.inertia## inertia cum cum(%)
## Ax1 7.57389986 7.57390 63.11583
## Ax2 2.88991089 10.46381 87.19842
## Ax3 0.60177036 11.06558 92.21318
## Ax4 0.38830444 11.45389 95.44905
## Ax5 0.24642988 11.70032 97.50263
## Ax6 0.15677898 11.85709 98.80912
## Ax7 0.07184318 11.92894 99.40781
## Ax8 0.04578521 11.97472 99.78936
## Ax9 0.02527719 12.00000 100.00000La tabla de correlacion corresponde a la primera tabla, de la podemos ver como DA esta relacionado fuerte y en el mismo sentido con el color, por otro lado, el aroma esta relacionado pero en sentido opuesto al color, es decir, entre mas claro el cafe mas suave el aroma, las demas variables se interpretan de manera similar.
La segunda tabla me permite identifacar la inercia, en la cual nos damos cuenta cuantas variables escogemos de la descomposicion espectral, para este caso con las 2 primeras componentes denemos un 87,2% de informacion de toda la base de datos, si escogemos 3, ya tendriamos un 92,2% de informacion de toda la base de los cafes.
Ahora se realizar el circulo de correlaciones, que tiene una interpretacion parecida a la explicada anteriormente, sin embargo este nos permite visualizar la informacion.
###### Circulo de Correlaciones ######
s.corcircle(acp_cafe$co)
Note que al igual que la correlaciones, si estan muy juntas significa que estan muy relacionas en el mismo sentido, por ejemplo, el cuerpo, amargo y aroma estan muy cerca uno del otro, por lo tanto un cafe muy amargo tendra mucho cuerpo y un aroma bastante fuerte. Por otro lado, si se encuentran opuestos en el circulo de correlaciones, significan que estan muy relacionados pero en sentido opuesto, por ejemplo, el amargo y el color estan muy relacionados en sentido opuesto, es decir, si el cafe es muy oscuro entonces es mas amargo.
Finalmente veamos el grafico de los individuos.
###### Grafico de los Individuos ######
s.label(acp_cafe$li)
tipo=c("Cl","M","C","M","C","O","M","C","M","C")
tipo=as.factor(tipo)
s.class(acp_cafe$li,tipo)
Note que los excelsos tanto claro como oscuro, cuentan con las mejores caracteristicas, por otro laso los que tiene 20%- 40% cevada tienen el mejor color.
PUNTO 4: Clasificacion usando FactoClass
Sobre los mismos datos que trabajamos en las islas galápagos realizaremos un análisis de conglomerados que permita identificar las islas que son mas similares y que características tienen dichas islas. Por lo que realizaremos todo el procedimiento, que va desde identificar el valor mas apropiado de k (total de clúster) hasta la caracterización de los grupos por medio de pruebas de diferencia de medias t.
Realizamos tres grupos los cuales nos permitieron clasificar las diferentes islas, y todo esto segun las caracterizticas similares con las que estas contaban.
## The number of retained axes for factorial analysis is 3
##
## The number of axes for clustering is 3## Look the histogram of 25 indexes
## Partition in 3 clustersNote que con factoclass realizamos una clasificacion de los datos.
cluster$cluster## Baltra Bartolome Caldwell Champion Coamano Daphne.Major
## 1 1 1 1 1 1
## Daphne.Minor Darwin Eden Enderby Espanola Fernandina
## 1 2 1 1 1 3
## Gardner1 Gardner2 Genovesa Isabela Marchena Onslow
## 1 1 2 3 2 1
## Pinta Pinzon Las.Plazas Rabida SanCristobal SanSalvador
## 2 1 1 1 3 3
## SantaCruz SantaFe SantaMaria Seymour Tortuga Wolf
## 3 1 3 1 1 2
## Levels: 1 2 3cluster$carac.cont## class: 1
## Test.Value Class.Mean Frequency Global.Mean
## Area -2.120 7.220 19 261.709
## Scruz -2.644 31.984 19 56.977
## Nearest -3.098 3.916 19 10.060
## Species -3.221 33.947 19 85.233
## Endemics -3.450 13.000 19 26.100
## Elevation -3.501 163.000 19 368.033
## ------------------------------------------------------------
## class: 2
## Test.Value Class.Mean Frequency Global.Mean
## Nearest 4.238 34.76 5 10.060
## Scruz 4.016 168.52 5 56.977
## ------------------------------------------------------------
## class: 3
## Test.Value Class.Mean Frequency Global.Mean
## Endemics 4.699 73.000 6 26.100
## Species 4.677 281.000 6 85.233
## Elevation 4.459 1054.500 6 368.033
## Area 3.133 1250.417 6 261.709Note que las islas del primer grupo se esncuentran por debajo de la media global en casi todas las variables analizadas, sin embargo las del grupo dos y tres estan por encima o igual de la media global en todas sus variables, lo que nos permite concluir que las islas del grupo uno tienen menos area que las demas, ademas menos numero de especies de plantas, las del grupo dos estan en la media y las del grupo tres si cuentan con mas area, ademas tienen mayor numero de especies.