LAB2
Paula Javierre
2022-10-18
Estadística descriptiva y gráficos con R y RCommander
Ejercicio 1 (ejercicio guiado)
En este ejercicio construimos las variables en nuestro script, recodificamos variables no numéricas y buscamos resúmenes estadísticos y tablas de frecuencias de nuestras variables y datos:
# Definición de los variables en vectores
id=c(100,110,120,130,140,150,160,170,180,190)
edad=c(18,19,NA,18,24,27,22,25,22,15)
sexo=c(2,1,2,2,1,2,1,1,2,1)
peso=c(65,58,56,61,84,99,50,64,87,87)
altura=c(161,170,174,165,150,171,181,170,184,190)
niv_col=c(1,2,3,1,3,2,3,1,2,3)
# Definimos los factores con sus etiquetas (recodificamos variables)
sexo=factor(sexo, levels=c(1,2), labels=c("Hombre","Mujer"))
niv_col=factor(niv_col, levels=c(1,2,3), labels=c("Colesterol Bajo", "Colesterol Normal","Colesterol Alto"))
# Ordenamos las variables
sort(edad)## [1] 15 18 18 19 22 22 24 25 27# Diferentes formas de buscar estadísticos
summary(edad)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 15.00 18.00 22.00 21.11 24.00 27.00 1fivenum(peso)## [1] 50.0 58.0 64.5 87.0 99.0mean(peso)## [1] 71.1var(altura)## [1] 134.9333cor(altura,peso)## [1] 0.07103622sd(altura)## [1] 11.61608Ejercicio 2 (ejercicio guiado)
Se quiere estudiar la posible asociación entre el hecho de que una gestante fume durante el embarazo y que el bebé presente bajo peso al nacer. Se realiza un estudio de seguimiento de 2.000 embarazadas. Los resultados de este estudio se muestran a continuación: http://www.fisterra.com/mbe/investiga/chi/chi.asp.
# Creamos nuestra tabla de contingencia con los datos del estudio
Si<-c(43,105)
No<-c(207,1645)
tabla1<-data.frame(Si,No) # Categorías de Bajo Peso al Nacer
rownames(tabla1)<-c("Fuma","No fuma") # Categorías de Fumar en el Embarazo
tabla1## Si No
## Fuma 43 207
## No fuma 105 1645#Hacemos el test Chi-cuadrado para ver la relación o independencia de las variables
chisq.test(tabla1)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabla1
## X-squared = 38.427, df = 1, p-value = 5.685e-10Al tener un p-valor < 0.05, podemos llegar a la conclusión de que hay una asociación significativa entre las variables.
Ejercicio 3 (ejercicio para practicar)
¡Ahora vosotros!: Escogemos dos conjuntos de datos incorporados en los paquetes Datasets y MASS de R. Los data frame que queremos son Orange. <(https://stat.ethz.ch/R-manual/R-atched/library/datasets/html/Orange.html)> e Iris del que ya hemos hablado anteriormente.
Deberemos:
1. Buscar un resumen de las variables de cada dataset.
# Cargamos los paquetes
data(package="MASS")
data(package="datasets")
data("Orange")
data("iris")
# Extraemos los datos más importantes de cada paquete
# options(max.print = 20)
head(Orange)## Tree age circumference
## 1 1 118 30
## 2 1 484 58
## 3 1 664 87
## 4 1 1004 115
## 5 1 1231 120
## 6 1 1372 142tail(Orange)## Tree age circumference
## 30 5 484 49
## 31 5 664 81
## 32 5 1004 125
## 33 5 1231 142
## 34 5 1372 174
## 35 5 1582 177names(Orange)## [1] "Tree" "age" "circumference"dim(Orange)## [1] 35 3summary(Orange)## Tree age circumference
## 3:7 Min. : 118.0 Min. : 30.0
## 1:7 1st Qu.: 484.0 1st Qu.: 65.5
## 5:7 Median :1004.0 Median :115.0
## 2:7 Mean : 922.1 Mean :115.9
## 4:7 3rd Qu.:1372.0 3rd Qu.:161.5
## Max. :1582.0 Max. :214.0head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosatail(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 145 6.7 3.3 5.7 2.5 virginica
## 146 6.7 3.0 5.2 2.3 virginica
## 147 6.3 2.5 5.0 1.9 virginica
## 148 6.5 3.0 5.2 2.0 virginica
## 149 6.2 3.4 5.4 2.3 virginica
## 150 5.9 3.0 5.1 1.8 virginicanames(iris)## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"dim(iris)## [1] 150 5summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## - Generar una tabla de frecuencias absolutas y una de relativas con cualquier variable del dataset.
# Miramos la frecuencia absoluta y relativa entre la longitud del pétalo y la especie
table(iris$Petal.Length,iris$Species) # Absoluta##
## setosa versicolor virginica
## 1 1 0 0
## 1.1 1 0 0
## 1.2 2 0 0
## 1.3 7 0 0
## 1.4 13 0 0
## 1.5 13 0 0
## 1.6 7 0 0
## 1.7 4 0 0
## 1.9 2 0 0
## 3 0 1 0
## 3.3 0 2 0
## 3.5 0 2 0
## 3.6 0 1 0
## 3.7 0 1 0
## 3.8 0 1 0
## 3.9 0 3 0
## 4 0 5 0
## 4.1 0 3 0
## 4.2 0 4 0
## 4.3 0 2 0
## 4.4 0 4 0
## 4.5 0 7 1
## 4.6 0 3 0
## 4.7 0 5 0
## 4.8 0 2 2
## 4.9 0 2 3
## 5 0 1 3
## 5.1 0 1 7
## 5.2 0 0 2
## 5.3 0 0 2
## 5.4 0 0 2
## 5.5 0 0 3
## 5.6 0 0 6
## 5.7 0 0 3
## 5.8 0 0 3
## 5.9 0 0 2
## 6 0 0 2
## 6.1 0 0 3
## 6.3 0 0 1
## 6.4 0 0 1
## 6.6 0 0 1
## 6.7 0 0 2
## 6.9 0 0 1prop.table(table(iris$Petal.Length,iris$Species)) # Relativa##
## setosa versicolor virginica
## 1 0.006666667 0.000000000 0.000000000
## 1.1 0.006666667 0.000000000 0.000000000
## 1.2 0.013333333 0.000000000 0.000000000
## 1.3 0.046666667 0.000000000 0.000000000
## 1.4 0.086666667 0.000000000 0.000000000
## 1.5 0.086666667 0.000000000 0.000000000
## 1.6 0.046666667 0.000000000 0.000000000
## 1.7 0.026666667 0.000000000 0.000000000
## 1.9 0.013333333 0.000000000 0.000000000
## 3 0.000000000 0.006666667 0.000000000
## 3.3 0.000000000 0.013333333 0.000000000
## 3.5 0.000000000 0.013333333 0.000000000
## 3.6 0.000000000 0.006666667 0.000000000
## 3.7 0.000000000 0.006666667 0.000000000
## 3.8 0.000000000 0.006666667 0.000000000
## 3.9 0.000000000 0.020000000 0.000000000
## 4 0.000000000 0.033333333 0.000000000
## 4.1 0.000000000 0.020000000 0.000000000
## 4.2 0.000000000 0.026666667 0.000000000
## 4.3 0.000000000 0.013333333 0.000000000
## 4.4 0.000000000 0.026666667 0.000000000
## 4.5 0.000000000 0.046666667 0.006666667
## 4.6 0.000000000 0.020000000 0.000000000
## 4.7 0.000000000 0.033333333 0.000000000
## 4.8 0.000000000 0.013333333 0.013333333
## 4.9 0.000000000 0.013333333 0.020000000
## 5 0.000000000 0.006666667 0.020000000
## 5.1 0.000000000 0.006666667 0.046666667
## 5.2 0.000000000 0.000000000 0.013333333
## 5.3 0.000000000 0.000000000 0.013333333
## 5.4 0.000000000 0.000000000 0.013333333
## 5.5 0.000000000 0.000000000 0.020000000
## 5.6 0.000000000 0.000000000 0.040000000
## 5.7 0.000000000 0.000000000 0.020000000
## 5.8 0.000000000 0.000000000 0.020000000
## 5.9 0.000000000 0.000000000 0.013333333
## 6 0.000000000 0.000000000 0.013333333
## 6.1 0.000000000 0.000000000 0.020000000
## 6.3 0.000000000 0.000000000 0.006666667
## 6.4 0.000000000 0.000000000 0.006666667
## 6.6 0.000000000 0.000000000 0.006666667
## 6.7 0.000000000 0.000000000 0.013333333
## 6.9 0.000000000 0.000000000 0.006666667- Generar una tabla de frecuencias absolutas con cada una de las variables del conjunto de datos Orange. ¿Todas las tablas generadas tienen sentido para vosotros?
# Miramos la frecuencia absoluta entre la edad y la circumferencia de dataset Orange
table(Orange$age,Orange$circumference) # Absoluta##
## 30 32 33 49 51 58 62 69 75 81 87 108 111 112 115 120 125 139 140 142 145
## 118 3 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 484 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 664 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0 0 0 0 0 0 0
## 1004 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0
## 1231 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0
## 1372 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0
## 1582 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
##
## 156 167 172 174 177 179 203 209 214
## 118 0 0 0 0 0 0 0 0 0
## 484 0 0 0 0 0 0 0 0 0
## 664 0 0 0 0 0 0 0 0 0
## 1004 1 1 0 0 0 0 0 0 0
## 1231 0 0 1 0 0 1 0 0 0
## 1372 0 0 0 1 0 0 1 1 0
## 1582 0 0 0 0 1 0 1 0 1# Miramos la frecuencia absoluta entre el árbol y la circumferencia de dataset Orange
table(Orange$Tree,Orange$circumference) # Absoluta##
## 30 32 33 49 51 58 62 69 75 81 87 108 111 112 115 120 125 139 140 142 145
## 3 1 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 1 0 0
## 1 1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 1
## 5 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0
## 2 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0
## 4 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0
##
## 156 167 172 174 177 179 203 209 214
## 3 0 0 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0
## 5 0 0 0 1 1 0 0 0 0
## 2 1 0 1 0 0 0 2 0 0
## 4 0 1 0 0 0 1 0 1 1# Miramos la frecuencia absoluta entre el árbol y la edad de dataset Orange
table(Orange$Tree,Orange$age) # Absoluta##
## 118 484 664 1004 1231 1372 1582
## 3 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1Comparando las 3 tablas, vemos que la que más sentido cobra es la que compara el árbol con la edad de las naranjas. Las otras dos tablas también muestran relaciones interesantes pero es difícil interpretar los resultados y encontrar una tendencia.
Ejercicio 4 (ejercicio guiado)
En este ejercicio realizaremos paso a paso la creación de gráficos y modificación de sus opciones para que veáis cómo se pueden generar:
library (MASS)
data("iris")
#Generamos los gráficos que hemos resumido en la tabla anterior
stem(iris$Sepal.Length)##
## The decimal point is 1 digit(s) to the left of the |
##
## 42 | 0
## 44 | 0000
## 46 | 000000
## 48 | 00000000000
## 50 | 0000000000000000000
## 52 | 00000
## 54 | 0000000000000
## 56 | 00000000000000
## 58 | 0000000000
## 60 | 000000000000
## 62 | 0000000000000
## 64 | 000000000000
## 66 | 0000000000
## 68 | 0000000
## 70 | 00
## 72 | 0000
## 74 | 0
## 76 | 00000
## 78 | 0# The stem and leaf plot is very useful for displaying the shape and relative density of data, hence giving the reader or customer a quick overview of the kind of distribution. The stem() command extracts the numeric data and splits them into two parts namely, the stem and leaf.
boxplot(iris$Petal.Width)
hist(iris$Sepal.Length)
plot(iris$Sepal.Width)
#Podemos cambiar opciones de los gráficos para hacer que tengan un mejor aspecto #Modificamos la escala y la medida del espacio gráfico del gráfico de tallo y hojas
stem(iris$Sepal.Length, scale =1, width =80)##
## The decimal point is 1 digit(s) to the left of the |
##
## 42 | 0
## 44 | 0000
## 46 | 000000
## 48 | 00000000000
## 50 | 0000000000000000000
## 52 | 00000
## 54 | 0000000000000
## 56 | 00000000000000
## 58 | 0000000000
## 60 | 000000000000
## 62 | 0000000000000
## 64 | 000000000000
## 66 | 0000000000
## 68 | 0000000
## 70 | 00
## 72 | 0000
## 74 | 0
## 76 | 00000
## 78 | 0#Añadimos color rojo al diagrama de caja y un título
boxplot(iris$Petal.Width, col="tomato", title="Diagrama de caja")
#Añadimos los intervalos que creamos pertinentes y le damos color azul al histograma
hist(iris$Sepal.Length, breaks=20, col="cadetblue3")
#Añadimos un título al eje de las "x", al eje de las "y" y el color rojo oscuro
plot (iris$Sepal.Width, col="dark red", xlab="SepalWidth", ylab="Y")
Ejercicio 5 (ejercicio de práctica)
¡Ahora vosotros!: Seguid los pasos y solucionad el ejercicio.
- Incorporad el dataset VaDeaths de R: https://www.rdocumentation.org/packages/datasets/versions/3.6.1/topics/VADeaths) y generad un gráfico barplot.
data("VADeaths") # Death rates per 1000 in Virginia in 1940.
VADeaths## Rural Male Rural Female Urban Male Urban Female
## 50-54 11.7 8.7 15.4 8.4
## 55-59 18.1 11.7 24.3 13.6
## 60-64 26.9 20.3 37.0 19.3
## 65-69 41.0 30.9 54.6 35.1
## 70-74 66.0 54.3 71.1 50.0# The death rates are measured per 1000 population per year. They are cross-classified by age group (rows) and population group (columns). The age groups are: 50--54, 55--59, 60--64, 65--69, 70--74 and the population groups are Rural/Male, Rural/Female, Urban/Male and Urban/Female.
# Cargamos una librería para poner color a los gráficos
library("RColorBrewer")
# display.brewer.all()
barplot(VADeaths,xlab="Population group",ylab="Number of deaths",main="Death rates in Virginia in 1940",col=brewer.pal(n = 5, name = "RdYlBu"),legend.text=TRUE)
- Usad la función pairs sobre el conjunto de datos iris.
Para este tipo de gráficos, se deben interpretar cada uno de ellos como la intersección entre la variable columna y fila. Sin embargo, al ser tan pequeños y con tantos datos, es complicado sacar conclusiones.
# ?pairs() --> A matrix of scatterplots is produced.
# Se puede introducir solamente iris o acotar las columnas que queremos incluir en el scatter plot
pairs(iris[1:5], main = "Scatterplot de iris",
pch = 21, bg = c("red", "yellow", "blue")[unclass(iris$Species)])
- Generad unos datos inventados y cread un boxplot con
ellos.
Utilizaremos la función intengrada runif. Esta nos permite generar un número determinado de muestra aleatorias entre los dos valores que especifiquemos
# Generamos datos aleatorios referentes a las edades de una población. Añadiendo floor nos aseguramos de tener número enteros.
randomdata <- floor(runif(100, 0, 90)) # Para crear valores random
randomdata## [1] 9 53 62 76 11 81 80 54 71 55 64 30 6 31 49 78 38 31 0 58 89 23 68 46 65
## [26] 81 28 33 19 34 30 33 56 31 80 12 17 83 6 48 0 0 51 46 28 39 47 39 62 3
## [51] 78 79 45 66 60 12 68 70 88 31 62 14 65 84 71 10 80 51 36 72 1 43 75 24 73
## [76] 38 42 0 54 4 71 57 1 37 45 41 57 39 60 76 5 61 85 40 57 81 78 44 31 34boxplot(randomdata,col=brewer.pal(n = 1, name = "Pastel2") )## Warning in brewer.pal(n = 1, name = "Pastel2"): minimal value for n is 3, returning requested palette with 3 different levels
- Dibujad una parábola y=x^2 con valores que van de x entre -10 y 10.
randomdata2 <- seq(-10,10) # Crea vector con valores de -10 a 10 con un step de 1
randomdata2## [1] -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8
## [20] 9 10parabola <- randomdata2^2
parabola## [1] 100 81 64 49 36 25 16 9 4 1 0 1 4 9 16 25 36 49 64
## [20] 81 100plot(randomdata2,parabola,xlab="x",ylab="x^2",type="l",main="Parábola")
Ejercicio 6 (ejercicio de práctica)
¡Ahora vosotros!: Cargad (o instalad primero y luego cargad) el paquete UsingR (https://cran.r-project.org/web/packages/UsingR/UsingR.pdf). El conjunto de datos brightness (incorporado en este paquete) contiene información sobre el brillo de 963 estrellas.
- Representad estos datos mediante un histograma y un gráfico de densidad superpuesto. Combinad los dos gráficos.
# Instalamos el paquete UsingR si no está instalado
# install.packages("UsingR")
# Cargamos el paquete
library(UsingR)
data(brightness)
dbrig<-brightness
hist(dbrig,prob=TRUE,main="Densidad de brillo de 963 estrellas",col="bisque",ylim=c(0,0.5))
lines(density(dbrig),col="coral")
- Representad gráficamente estos datos mediante un diagrama de caja (boxplot).
boxplot(dbrig,data=dbrig,main="Brillo de 963 estrellas", xlab="% estrellas", ylab="Valores", col="darkseagreen3")
- ¿Diríais que los datos presentan «outliers»?
Un outlier es un valor o una observación que está distante de otras observaciones, es decir, un punto de datos que difiere significativamente de otros puntos de datos. Estos puntos son los que vemos en el gráfico del apartado anterior.Podemos concluir que sí, existen outliers en nuestros datos.
Ejercicio 7 (ejercicio de práctica)
¡Ahora vosotros!: Representad gráficamente los datos contenidos en los conjuntos de datos (incorporados en UsingR) bumpers, firstchi, math con un histograma y/o un boxplot usando ggplot2.
# Cargamos los paquetes
data(bumpers)
data(firstchi)
data(math)
library(ggplot2)
?bumpers # Price in dollars to repair a bumper
?firstchi # Age of mother at birth of first child
?math # Standardized math scores
# Para los datos bumpers
qplot(bumpers, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Price in dollars to repair a bumper",xlab="Number of bumpers",ylab="Dollars")
boxplot(bumpers, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Price in dollars to repair a bumper",xlab="Number of bumpers",ylab="Dollars")
# Para los datos firstchi
qplot(firstchi, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Age of mother at birth of first child",xlab="Number of mothers",ylab="Age")
boxplot(firstchi, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Age of mother at birth of first child",xlab="Number of mothers",ylab="Age")
# Para los datos math
qplot(math, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Standardized math scores",xlab="Repetitions",ylab="Score")
boxplot(math, geom ='histogram', fill=I("lightskyblue"), col=I("cornflowerblue"), main="Standardized math scores",xlab="Repetitions",ylab="Score")
Ejemplo de análisis de regresión lineal simple paso a paso
Leemos el fichero Agro.csv* y generamos un dataset con el nombre Agro. El fichero lo encontraréis adjunto en el LAB. Se pretende ver la relación que hay entre la cantidad de agua aplicada (en m3) y el correspondiente rendimiento de una cosecha de fresas (en toneladas por hectárea).
#P1: Cargamos los datos, observamos el dataset y resumimos nuestras variables.
library(readxl)
Agro <-read_excel("G:/Mi unidad/SW/LAB2_Data/LAB2_Data/Agro.xlsx")
View(Agro)
# Nombre de las variables
names(Agro)## [1] "Cosecha" "Agua"# Resumen de las variables
summary(Agro)## Cosecha Agua
## Min. : 8 Min. : 4.10
## 1st Qu.:30 1st Qu.: 6.85
## Median :52 Median : 9.69
## Mean :52 Mean : 9.69
## 3rd Qu.:74 3rd Qu.:12.49
## Max. :96 Max. :15.30#P2: Representamos las variables con un gráfico de dispersión usando la función plot(). En nuestro caso, la variable Agua será la variable "x" y la variable Cosecha, la "y".
plot(Agro$Cosecha,Agro$Agua)
#P3: Como vemos en este gráfico, existe una gran relación lineal entre la cantidad de agua que usamos y las toneladas de fresas que recogemos en la cosecha. En este punto generamos nuestro modelo lineal y estimamos los parámetros que tendrá.
modelo<-lm(Agro$Cosecha~Agro$Agua, data=Agro)
summary(modelo)##
## Call:
## lm(formula = Agro$Cosecha ~ Agro$Agua, data = Agro)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.314040 -0.000323 0.000430 0.001022 0.313179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -23.99796 0.13259 -181.0 <2e-16 ***
## Agro$Agua 7.84293 0.01286 609.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1569 on 10 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.719e+05 on 1 and 10 DF, p-value: < 2.2e-16#P4: Para ajustar la recta de regresión debemos mirar el valor estimado de los parámetros en el resumen anterior. Nos fijamos en el coeficiente de determinación (R2) o "R-squared" para comprobar cómo de bueno es el ajuste. Es 1, por lo que es completamente un modelo lineal.
#P5: Miramos la correlación de Pearson con la matriz de coeficientes de correlación y con el test de Pearson.
cor(Agro)## Cosecha Agua
## Cosecha 1.0000000 0.9999866
## Agua 0.9999866 1.0000000cor.test(Agro$Cosecha, Agro$Agua, method="pearson")##
## Pearson's product-moment correlation
##
## data: Agro$Cosecha and Agro$Agua
## t = 609.87, df = 10, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9999503 0.9999964
## sample estimates:
## cor
## 0.9999866#P6: Miramos la normalidad, es decir, si los errores estandarizados siguen una distribución normal con 3 gráficos.
residuos<-rstandard(modelo) # residuos estándares del modelo ajustado (completo)
par(mfrow=c(1,3)) # Para mostrar los 3 gráficos en 1, creamos la matriz necesaria
hist(residuos) # histograma de los residuos estandarizados
boxplot(residuos) # diagrama de cajas de los residuos estandarizados
qqnorm(residuos) # gráfico de cuantiles de los residuos estandarizados
qqline(residuos) # dibujamos la recta en el gráfico anterior
par(mfrow=c(1,1)) # devuelve la pantalla a su estado original sin divisione
#P7: Miramos si la varianza de los errores es constante con el gráfico siguiente:
plot(fitted.values(modelo),rstandard(modelo), xlab="Valores ajustados", ylab="Residuos estandarizados") #Generamos el gráfico
abline(h=0) # dibujamos la recta en el cero
#P8: Buscar valores atípico
plot(Agro$Agua,rstandard(modelo),xlab="Agua",ylab="Residuos estandarizados")
#P9: Visualizamos la regresión linea
plot(Agro$Agua, Agro$Cosecha ,xlab ="Cantidad de Agua (m3)", ylab ="Cosecha de fresas en toneladas" )
abline(modelo)
#P10: Utilizando regresión lineal, obtenemos el rendimiento que cabe esperar si la cantidad de agua aplicada es 5 m3. Podemos usar la función predict() o sustituir directamente en el modelo con parámetros estimados.
## (Intercept) -23.99796
## Agro$Agua 7.84293
predic=-23.99+7.84*(5)
predic## [1] 15.21Ejercicio 9: Ejercicio final R
A partir de unos datos bioclínicos o biosanitarios que escojáis y que importéis a R, explicad sus variables (un mínimo de 8) y:
- Realizad un resumen estadístico completo del dataset y explicad los resultados.
Utilizaremos el dataset Pima.te, el cual se encuentra dentro de la librería MASS y contiene valores de diabetes para una población determinada.
?Pima.te
- npreg: number of pregnancies.
- glu: plasma glucose concentration in an oral glucose tolerance
test.
- bp: diastolic blood pressure (mm Hg).
- skin: triceps skin fold thickness (mm).
- bmi: body mass index (weight in kg/(height in m)^2).
- ped: diabetes pedigree function.
- age: age in years.
- type: Yes or No, for diabetic according to WHO criteria.
library(MASS) # Cargamos la libreria MASS
data("Pima.te") # Cargamos el dataset que nos interesa
head(Pima.te)## npreg glu bp skin bmi ped age type
## 1 6 148 72 35 33.6 0.627 50 Yes
## 2 1 85 66 29 26.6 0.351 31 No
## 3 1 89 66 23 28.1 0.167 21 No
## 4 3 78 50 32 31.0 0.248 26 Yes
## 5 2 197 70 45 30.5 0.158 53 Yes
## 6 5 166 72 19 25.8 0.587 51 Yessummary(Pima.te)## npreg glu bp skin
## Min. : 0.000 Min. : 65.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.000 1st Qu.: 96.0 1st Qu.: 64.00 1st Qu.:22.00
## Median : 2.000 Median :112.0 Median : 72.00 Median :29.00
## Mean : 3.485 Mean :119.3 Mean : 71.65 Mean :29.16
## 3rd Qu.: 5.000 3rd Qu.:136.2 3rd Qu.: 80.00 3rd Qu.:36.00
## Max. :17.000 Max. :197.0 Max. :110.00 Max. :63.00
## bmi ped age type
## Min. :19.40 Min. :0.0850 Min. :21.00 No :223
## 1st Qu.:28.18 1st Qu.:0.2660 1st Qu.:23.00 Yes:109
## Median :32.90 Median :0.4400 Median :27.00
## Mean :33.24 Mean :0.5284 Mean :31.32
## 3rd Qu.:37.20 3rd Qu.:0.6793 3rd Qu.:37.00
## Max. :67.10 Max. :2.4200 Max. :81.00dim(Pima.te) # Nos aseguramos que tega mínimo 8 variables## [1] 332 8names(Pima.te)## [1] "npreg" "glu" "bp" "skin" "bmi" "ped" "age" "type"# Valores estadísticos
sapply(Pima.te, mean, na.rm=TRUE) # Los valores medios de cada variable. Los que no es aplicable, aparecen como NA## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA## npreg glu bp skin bmi ped
## 3.4849398 119.2590361 71.6536145 29.1626506 33.2397590 0.5283886
## age type
## 31.3162651 NAfivenum(Pima.te$age) # min, Q1, mean, Q3, max para las edades## [1] 21 23 27 37 81max(Pima.te$bmi) # BMI mínimo## [1] 67.1min(Pima.te$bmi) # BMI máximo## [1] 19.4var(Pima.te$bp) # Varianza en la presión ## [1] 163.8223Podemos utilizar también tablas para resumir datos relevantes. Por ejemplo:
table(Pima.te$age,Pima.te$npreg) # Relación edad con número de embarazos##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 15 17
## 21 9 11 5 1 1 0 0 0 0 0 0 0 0 0 0 0
## 22 10 12 11 6 2 0 0 0 0 0 0 0 0 0 0 0
## 23 5 7 4 1 1 0 1 0 0 0 0 0 0 0 0 0
## 24 5 7 7 2 2 0 0 0 0 0 0 0 0 0 0 0
## 25 7 7 10 5 1 1 0 0 0 0 0 0 0 0 0 0
## 26 4 4 1 4 2 0 1 0 0 0 0 0 0 0 0 0
## 27 2 2 2 3 1 2 0 0 0 0 0 0 0 0 0 0
## 28 0 6 1 5 1 0 3 0 0 0 0 0 0 0 0 0
## 29 1 4 3 0 1 1 0 1 0 1 0 0 0 0 0 0
## 30 0 1 1 2 1 3 0 0 0 0 0 0 0 0 0 0
## 31 1 2 0 2 1 0 0 0 1 0 0 0 0 0 0 0
## 32 1 0 0 0 0 2 2 1 0 0 0 0 0 0 0 0
## 33 0 2 1 0 3 0 1 0 0 0 0 0 0 0 0 0
## 34 0 0 1 0 1 0 1 1 0 1 0 0 0 0 0 0
## 35 1 0 0 1 0 2 1 0 0 0 0 0 0 0 0 0
## 36 0 2 0 2 1 0 0 2 2 0 0 0 0 0 0 0
## 37 0 1 0 0 2 1 1 1 1 0 0 0 0 0 0 0
## 38 1 1 0 1 0 1 0 0 0 1 2 0 0 0 0 0
## 39 0 0 0 1 0 0 1 1 2 0 0 0 0 1 0 0
## 40 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 0
## 41 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0
## 42 0 0 1 0 0 1 0 2 1 1 0 2 0 0 0 0
## 43 0 1 0 0 0 1 1 2 0 1 0 0 0 1 1 0
## 45 0 0 0 0 1 0 0 1 1 0 1 1 0 1 0 0
## 46 0 0 0 0 0 0 1 0 1 1 0 0 2 0 0 0
## 47 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 1
## 48 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0
## 49 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0
## 50 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0
## 51 0 0 0 0 0 1 0 2 0 0 0 1 0 0 0 0
## 52 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## 53 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0
## 54 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## 56 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0
## 57 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## 58 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 60 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0
## 61 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 63 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
## 65 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## 70 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## 81 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0# table(Pima.te$bmi,Pima.te$bp) # Relación del BMI con la presión arterial
table(Pima.te$glu,Pima.te$type) # Relación enttre concentración de glucosa y criterio para diabetes##
## No Yes
## 65 1 0
## 68 3 0
## 71 1 0
## 72 1 0
## 73 1 0
## 74 2 0
## 75 1 0
## 77 1 0
## 78 1 1
## 80 2 0
## 81 5 0
## 82 2 0
## 83 2 0
## 84 7 1
## 85 1 1
## 87 4 0
## 88 7 1
## 89 5 0
## 90 5 2
## 91 3 0
## 92 2 0
## 93 4 1
## 94 5 0
## 95 6 1
## 96 4 0
## 97 3 1
## 98 3 0
## 99 8 0
## 100 7 1
## 101 4 1
## 102 5 2
## 103 6 0
## 104 1 2
## 105 4 1
## 106 5 1
## 107 2 1
## 108 8 0
## 109 2 5
## 111 4 3
## 112 5 1
## 113 3 1
## 114 1 0
## 115 1 0
## 116 2 1
## 117 5 1
## 118 1 1
## 119 3 3
## 120 3 1
## 121 2 0
## 122 2 2
## 123 4 0
## 124 3 1
## 125 2 2
## 126 6 0
## 127 3 0
## 128 2 4
## 129 3 3
## 130 1 0
## 131 2 0
## 132 1 0
## 133 1 0
## 134 2 1
## 135 0 1
## 136 1 3
## 137 1 0
## 138 0 1
## 139 2 0
## 141 1 0
## 142 1 1
## 143 2 0
## 144 2 2
## 145 0 2
## 146 2 1
## 147 2 1
## 148 0 1
## 150 1 0
## 151 2 1
## 152 1 2
## 153 2 0
## 154 1 0
## 155 0 3
## 156 0 1
## 157 1 0
## 158 0 1
## 160 0 1
## 161 0 1
## 162 0 3
## 163 0 2
## 165 2 0
## 166 0 1
## 169 0 1
## 170 0 2
## 171 0 1
## 172 0 1
## 173 1 4
## 174 0 2
## 176 0 1
## 179 0 3
## 180 1 2
## 181 0 4
## 184 0 1
## 186 0 1
## 187 0 3
## 189 0 2
## 193 0 1
## 196 0 2
## 197 1 1prop.table(table(Pima.te$age,Pima.te$type)) # Frecuencia de diabéticos por edad##
## No Yes
## 21 0.078313253 0.003012048
## 22 0.111445783 0.012048193
## 23 0.048192771 0.009036145
## 24 0.060240964 0.009036145
## 25 0.063253012 0.030120482
## 26 0.030120482 0.018072289
## 27 0.024096386 0.012048193
## 28 0.033132530 0.015060241
## 29 0.018072289 0.018072289
## 30 0.015060241 0.009036145
## 31 0.006024096 0.015060241
## 32 0.012048193 0.006024096
## 33 0.012048193 0.009036145
## 34 0.012048193 0.003012048
## 35 0.006024096 0.009036145
## 36 0.012048193 0.015060241
## 37 0.012048193 0.009036145
## 38 0.009036145 0.012048193
## 39 0.015060241 0.003012048
## 40 0.006024096 0.006024096
## 41 0.003012048 0.006024096
## 42 0.018072289 0.006024096
## 43 0.006024096 0.018072289
## 45 0.009036145 0.009036145
## 46 0.006024096 0.009036145
## 47 0.003012048 0.006024096
## 48 0.006024096 0.003012048
## 49 0.006024096 0.003012048
## 50 0.000000000 0.009036145
## 51 0.003012048 0.009036145
## 52 0.000000000 0.003012048
## 53 0.000000000 0.009036145
## 54 0.000000000 0.003012048
## 56 0.003012048 0.006024096
## 57 0.000000000 0.003012048
## 58 0.003012048 0.000000000
## 60 0.006024096 0.000000000
## 61 0.003012048 0.000000000
## 63 0.006024096 0.000000000
## 65 0.003012048 0.000000000
## 70 0.000000000 0.003012048
## 81 0.003012048 0.000000000- Realizad 5 gráficos con las variables, explicad su significado y guardadlos como imágenes (jpeg o bmp).
Primero comentaremos los gráficos uno a uno:
En el primero podemos apreciar la frecuencia de embarazos en la
población de mujeres. Vemos que la gran mayoría de mujeres tienen 5 o
menos embarazos, pero hay individuos que tienen hasta 16.
hist(Pima.te$npreg,main="Number of pregnancies",xlab="Women",ylab="Frequency",ylim=c(0,200),xlim=c(0,20),col="darkorchid1")
Otro parámetro que nos interesa tener siempre en cuenta cuando analizamos poblaciones, es la edad de la muestra. Con un boxplot podemos obtener los valores mínimos, máximos, cuartiles y la mediana. Podemos ver que la edad mediana es aproximadamente 25. También es importante remarcar que tenemos outliers o valores atípicos, los cuales son representados por los círculos y nos informan de los valores que superan la varianza esperada.
boxplot(Pima.te$age,col=brewer.pal(n = 1, name = "Pastel2"),main="Age of pima Indian women",ylab="Age (years)")
Los diagramas de cajas también nos permiten comparar dos grupos de estudio diferentes para una variable común. De dicha forma podemos ver que en nuestra población, las mujeres diabéticas tienen un índice de masa corporal superior a las no diabéticas.
plot(Pima.te$type,Pima.te$bmi,col=brewer.pal(n = 2, name = "BuPu"),main="BMI for diabetic and non diabetic women",ylab="BMI (weight in kg/(height in m)^2) ",names=c("Non diabetic","Diabetic"),xlab="WHO criteria")
Un gráfico simple y sencillo es el siguiente. En él vemos la relación entre el índice de masa corporal y el grosor de la piel del tríceps. Al graficar los puntos vemos una tendencia lineal. Para poder compararlo, se ha dibujado encima la línea y=x, la cual nos permitiría poder asumir o al menos valorar que estas dos variables siguen una regresión lineal: a mayor BMI, mayor grosor de piel. Lo analizaremos con más detalle en el siguiente apartado.
plot(Pima.te$bmi,Pima.te$skin,col="blue2",main="Relation between BMI and triceps skin thicnkess",ylab="Triceps skin thickness (mm)",xlab="BMI (weight in kg/(height in m)^2",xlim=c(0,70),ylim=c(0,70))
x=seq(0,100)
lines(x,col="red")
También se ha querido revisar si las diabéticas tienen un mayor grosor de piel en el tríceps. Los resultados son esperados, ya que anteriormente hemos visto que las diabéticas tienen mayor BMI, y que las mujeres con mayor BMI son diabéticas.
plot(Pima.te$type,Pima.te$skin,col=brewer.pal(n = 2, name = "Spectral"),main="Skin thicnkess for diabetic and non diabetic women",ylab="Triceps skin fold thickness (mm)",names=c("Non diabetic","Diabetic"),xlab="WHO criteria")
# Para guardar las imágenes
jpeg(file="Number of pregnancies.jpeg")
hist(Pima.te$npreg,main="Number of pregnancies",xlab="Women",ylab="Frequency",ylim=c(0,200),xlim=c(0,20),col="darkorchid1")
dev.off()## png
## 2bmp(file="Age of pima Indian women.bmp")
boxplot(Pima.te$age,col=brewer.pal(n = 1, name = "Pastel2"),main="Age of pima Indian women",ylab="Age (years)")
dev.off()## png
## 2jpeg(file="BMI for diabetic and non diabetic women.jpeg")
plot(Pima.te$type,Pima.te$bmi,col=brewer.pal(n = 2, name = "BuPu"),main="BMI for diabetic and non diabetic women",ylab="BMI (weight in kg/(height in m)^2) ",names=c("Non diabetic","Diabetic"),xlab="WHO criteria")
dev.off()## png
## 2jpeg(file="BMI and triceps skin.bmp")
plot(Pima.te$bmi,Pima.te$skin,col="blue2",main="Relation between BMI and triceps skin thicnkess",ylab="Triceps skin thickness (mm)",xlab="BMI (weight in kg/(height in m)^2",xlim=c(0,70),ylim=c(0,70))
x=seq(0,100)
lines(x,col="red")
dev.off()## png
## 2jpeg(file="Skin thicnkess.jpeg")
plot(Pima.te$type,Pima.te$skin,col=brewer.pal(n = 2, name = "Spectral"),main="Skin thicnkess for diabetic and non diabetic women",ylab="Triceps skin fold thickness (mm)",names=c("Non diabetic","Diabetic"),xlab="WHO criteria")
dev.off()## png
## 2- Generad una regresión lineal entre 2 de sus variables paso a paso y comentad los resultados obtenidos.
Después de los resultados obtenidos en el apartado anterior, vamos a analizar la posible regresión lineal entre el grosor de la piel y el índice de masa corporal (BMI).
#P1: Cargamos los datos, observamos el dataset y resumimos nuestras variables.
# (Ya se ha hecho en el apartado anterior).
library(MASS) # Cargamos la libreria MASS
data("Pima.te") # Cargamos el dataset que nos interesa
View(Pima.te)
summary(Pima.te)## npreg glu bp skin
## Min. : 0.000 Min. : 65.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.000 1st Qu.: 96.0 1st Qu.: 64.00 1st Qu.:22.00
## Median : 2.000 Median :112.0 Median : 72.00 Median :29.00
## Mean : 3.485 Mean :119.3 Mean : 71.65 Mean :29.16
## 3rd Qu.: 5.000 3rd Qu.:136.2 3rd Qu.: 80.00 3rd Qu.:36.00
## Max. :17.000 Max. :197.0 Max. :110.00 Max. :63.00
## bmi ped age type
## Min. :19.40 Min. :0.0850 Min. :21.00 No :223
## 1st Qu.:28.18 1st Qu.:0.2660 1st Qu.:23.00 Yes:109
## Median :32.90 Median :0.4400 Median :27.00
## Mean :33.24 Mean :0.5284 Mean :31.32
## 3rd Qu.:37.20 3rd Qu.:0.6793 3rd Qu.:37.00
## Max. :67.10 Max. :2.4200 Max. :81.00names(Pima.te)## [1] "npreg" "glu" "bp" "skin" "bmi" "ped" "age" "type"#P2: Representamos las variables con un gráfico de dispersión usando la función plot(). En nuestro caso, la variable skin será la variable "x" y la variable BMI, la "y".
plot(Pima.te$skin,Pima.te$bmi)
#P3: En este punto generamos nuestro modelo lineal y estimamos los parámetros que tendrá.
mod<-lm(Pima.te$bmi~Pima.te$skin, data=Pima.te)
summary(mod)##
## Call:
## lm(formula = Pima.te$bmi ~ Pima.te$skin, data = Pima.te)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.3118 -3.9398 -0.4137 3.3782 25.5762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.89158 0.95153 19.85 <2e-16 ***
## Pima.te$skin 0.49201 0.03095 15.90 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.489 on 330 degrees of freedom
## Multiple R-squared: 0.4337, Adjusted R-squared: 0.432
## F-statistic: 252.7 on 1 and 330 DF, p-value: < 2.2e-16#P4: Para ajustar la recta de regresión debemos mirar el valor estimado de los parámetros en el resumen anterior. El "R-squared" = 0.4337, lo cual nos indica que el ajuste no es bueno, no es completamente un modelo lineal. Habrá que valorar el modelo final.
#P5: Miramos la correlación para saber el grado de relación lineal. Lo hacemos con el test de Pearson.
cor.test(Pima.te$bmi, Pima.te$skin, method="pearson")##
## Pearson's product-moment correlation
##
## data: Pima.te$bmi and Pima.te$skin
## t = 15.897, df = 330, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5929347 0.7154654
## sample estimates:
## cor
## 0.6585428#P6: Miramos la normalidad, es decir, si los errores estandarizados siguen una distribución normal con 3 gráficos.
res<-rstandard(mod) # residuos estándares del modelo ajustado (completo)
par(mfrow=c(1,3)) # Para mostrar los 3 gráficos en 1, creamos la matriz necesaria
hist(res) # histograma de los residuos estandarizados
boxplot(res) # diagrama de cajas de los residuos estandarizados
qqnorm(res) # gráfico de cuantiles de los residuos estandarizados
qqline(res) # dibujamos la recta en el gráfico anterior
par(mfrow=c(1,1)) # devuelve la pantalla a su estado original sin divisiones
#P7: Miramos si la varianza de los errores es constante con el gráfico siguiente:
plot(fitted.values(mod),rstandard(mod), xlab="Valores ajustados", ylab="Residuos estandarizados") #Generamos el gráfico
abline(h=0) # dibujamos la recta en el cero
#P8: Buscamos valores atípico
plot(Pima.te$bmi,rstandard(mod),xlab="BMI",ylab="Residuos estandarizados")
#P9: Visualizamos la regresión ajustada a los puntos
plot(Pima.te$bmi, Pima.te$skin ,xlab ="Triceps skin fold thickness (mm)", ylab ="Body mass index (weight in kg/(height in m)^2)")
abline(mod)
#P10: Utilizando regresión lineal, obtenemos el BMI que cabe esperar si la skin thickness es 50 mm. Podemos usar la función predict() o sustituir directamente en el modelo con parámetros estimados.
## (Intercept) -18.89158
## Agro$Agua 0.49201
predic=18.90+0.49*(50)
predic## [1] 43.4Después de intentar realizar una regresión lineal entre las dos variables, podemos concluir que el modelo no se adapta perfectamente. Nuestro modelo explicaría en un 43.37% la variable dependiente. Es decir, el BMI solamente podría explicar o modelar ese porcentaje. Para poder ajustar el modelo, una teoría sería que este requiera de más variables a incluir.
Respecto al p-valor del modelo, diremos que es significativo, ya que es menor a p < 0.05. Este dato nos indica que el resultado de nuestro modelo es poco probable que sea aleatorio, sino que guarda relación.
Sabemos que el modelo de regresión lineal supone que los residuos siguen una distribución normal. Esta afirmación se cumple si revisamos los datos del summary del modelo y el gráfico cuantil-cuantil (Normal Q-Q). La mediana de los residuos es cercana a 0; y el 1er y 3er cuartil son prácticamente simétricos. En el gráfico Normal Q-Q comparamos la distribución de los residuos del modelo con la distribución normal.
Sin embargo, vemos en el histograma y en el boxplot que hay outliers y que la frecuencia de residuos es elevada; lo cual perjudica y perturba el modelo lineal. La varianza tampoco es simétrica, es decir, su media no es 0 y genera ruido en el modelo.
En conclusión, nuestro modelo de regresión lineal no está bien ajustado si solamente incluimos la variable citada. El valor simulado al final solamente lo podemos predecir con un 43,37% de seguridad. Habrá que valorar el introducir más variables al modelo y comprobar si se ajusta.