# The R Book
# capitulo 2
# 2.7.2 La funcion aggregate para sumario
# de estadisticas agrupadas
# Suponga que tenemos dos variables dependientes
# (y and z) y dos variables independientes (x and w)
# que nosotros podriamos desear usar para sumario de 2
# funciones como el promedio o la varianza de y and/or z.
# La funcion aggregate tiene un metodo en su formula
# el cual permite sumarios elegantes de 4 tipos:
# uno a uno aggregate(y ~ x, mean)
# one a muchas aggregate(y ~ x + w, mean)
# muchos a una aggregate(cbind(y,z) ~ x, mean)
# muchas a muchas aggregate(cbind(y,z) ~ x + w, mean)
# aqui se presenta un ejemplo usando una data frame
# con una variable continua Grow.rate y tres variables
# independientes (categoricas) (Water, Detergent y Daphnia)
# asignamos el archivo daphnia a la variable daphnia.data
setwd("C:\\Users\\Luis\\Documents\\therbook")
dapnia.data <- read.table("C:\\Users\\Luis\\Documents\\therbook\\daphnia.txt",
header = TRUE,stringsAsFactors = TRUE)
str(dapnia.data)
## 'data.frame': 72 obs. of 4 variables:
## $ Growth.rate: num 2.92 2.49 3.02 2.35 3.15 ...
## $ Water : Factor w/ 2 levels "Tyne","Wear": 1 1 1 1 1 1 1 1 1 1 ...
## $ Detergent : Factor w/ 4 levels "BrandA","BrandB",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ Daphnia : Factor w/ 3 levels "Clone1","Clone2",..: 1 1 1 2 2 2 3 3 3 1 ...
# uso de la funcion aggregate en el uso uno a uno
# para encontrar el promedio de Growth.rate en 2 las
# dos muestras de agua
aggregate(Growth.rate~Water,FUN= mean,data=dapnia.data)
## Water Growth.rate
## 1 Tyne 3.685862
## 2 Wear 4.017948
# Ejemplo de uno-a-muchas para ver la interaccion
# entre Water and Detergent
aggregate(Growth.rate~Water+Detergent,data=dapnia.data,
FUN=mean)
## Water Detergent Growth.rate
## 1 Tyne BrandA 3.661807
## 2 Wear BrandA 4.107857
## 3 Tyne BrandB 3.911116
## 4 Wear BrandB 4.108972
## 5 Tyne BrandC 3.814321
## 6 Wear BrandC 4.094704
## 7 Tyne BrandD 3.356203
## 8 Wear BrandD 3.760259
# 2.7.4 Sumario de informacion a partir de vectores por grupos
# La funcion vectorial tapply es una de las mas importantes
# y utiles . La ‘t’es por table la idea es aplicar la funcion
# para producir una tabla desde los valores del vector, basada
# en una variables de agrupamiento (usualmente el agrupamiento
# es por niveles del factor).
tapply(dapnia.data$Growth.rate,dapnia.data$Detergent, mean)
## BrandA BrandB BrandC BrandD
## 3.884832 4.010044 3.954512 3.558231
cancer <- read.table("C:\\Users\\Luis\\Documents\\therbook\\cancer.txt",
header = TRUE,stringsAsFactors = TRUE)
str(cancer)
## 'data.frame': 120 obs. of 3 variables:
## $ death : int 4 26 2 25 7 6 5 2 4 1 ...
## $ treatment: Factor w/ 4 levels "DrugA","DrugB",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ status : int 1 1 1 1 1 0 1 0 1 1 ...
aggregate(cancer$death~cancer$treatment,data = cancer,FUN=mean)
## cancer$treatment cancer$death
## 1 DrugA 9.633333
## 2 DrugB 8.500000
## 3 DrugC 6.266667
## 4 placebo 5.800000
germination.data <- read.table("C:\\Users\\Luis\\Documents\\therbook\\germination.txt",
header = TRUE,stringsAsFactors = TRUE)
str(germination.data)
## 'data.frame': 21 obs. of 4 variables:
## $ count : int 10 23 23 26 17 5 53 55 32 46 ...
## $ sample : int 39 62 81 51 39 6 74 72 51 79 ...
## $ Orobanche: Factor w/ 2 levels "a73","a75": 2 2 2 2 2 2 2 2 2 2 ...
## $ extract : Factor w/ 2 levels "bean","cucumber": 1 1 1 1 1 2 2 2 2 2 ...
aggregate(germination.data$sample~germination.data$Orobanche+
germination.data$extract,FUN=mean,
data = germination.data)
## germination.data$Orobanche germination.data$extract
## 1 a73 bean
## 2 a75 bean
## 3 a73 cucumber
## 4 a75 cucumber
## germination.data$sample
## 1 24.60000
## 2 54.40000
## 3 28.20000
## 4 49.16667
length(germination.data$sample)
## [1] 21
length(germination.data$Orobanche)
## [1] 21
tapply(germination.data$sample,germination.data$Orobanche,
FUN = mean)
## a73 a75
## 26.40000 51.54545
tapply(dapnia.data$Growth.rate, dapnia.data$Detergent,
FUN = mean)
## BrandA BrandB BrandC BrandD
## 3.884832 4.010044 3.954512 3.558231
# 2.7.5 Direcciones dentro de los vectores
# Con la funcion which es posible obtener informacion
# dentro del vector
y <- c(8,3,5,7,6,6,8,9,2,3,9,4,10,4,11)
which(y>5)
## [1] 1 4 5 6 7 8 11 13 15
# which devuelve un vector de indices, en este caso
# la posicion 1,4,5 etc
# Produce un vector mas cort que y
# porque los valores menor a 5 y 5 han
# quedado fuera
length(y)
## [1] 15
length(y[y>5])
## [1] 9
houses <- read.table("c:\\Users\\Luis\\Documents\\therbook\\houses.txt",header = TRUE)
names(houses)
## [1] "Location" "Price"
# Aplicamos tres funciones diferesntes
# al vector llamado Price:
ranks <- rank(houses$Price)
sorted <- sort(houses$Price)
ordered <- order(houses$Price)# ver pagina 49 the R-Book
view <- data.frame(houses$Price,ranks,sorted,ordered)
view
## houses.Price ranks sorted ordered
## 1 325 12.0 95 9
## 2 201 10.0 101 6
## 3 157 5.0 117 10
## 4 162 6.0 121 12
## 5 164 7.0 157 3
## 6 101 2.0 162 4
## 7 211 11.0 164 5
## 8 188 8.5 188 8
## 9 95 1.0 188 11
## 10 117 3.0 201 2
## 11 188 8.5 211 7
## 12 121 4.0 325 1
houses$Location[order(houses$Price)]#ver pagina 49
## [1] Reading Staines Winkfield Newbury Bracknell
## [6] Camberley Bagshot Maidenhead Warfield Sunninghill
## [11] Windsor Ascot
## 12 Levels: Ascot Bagshot Bracknell Camberley Maidenhead ... Winkfield
# matrices y arrays
y <- 1:24
dim(y) <- c(2,4,3)#la tercera dimension es 3
y
## , , 1
##
## [,1] [,2] [,3] [,4]
## [1,] 1 3 5 7
## [2,] 2 4 6 8
##
## , , 2
##
## [,1] [,2] [,3] [,4]
## [1,] 9 11 13 15
## [2,] 10 12 14 16
##
## , , 3
##
## [,1] [,2] [,3] [,4]
## [1,] 17 19 21 23
## [2,] 18 20 22 24
dim(y) <- c(3,2,4)#la tercera dimension es 4
y
## , , 1
##
## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6
##
## , , 2
##
## [,1] [,2]
## [1,] 7 10
## [2,] 8 11
## [3,] 9 12
##
## , , 3
##
## [,1] [,2]
## [1,] 13 16
## [2,] 14 17
## [3,] 15 18
##
## , , 4
##
## [,1] [,2]
## [1,] 19 22
## [2,] 20 23
## [3,] 21 24
vector <- c(1,2,3,4,4,3,2,1)
V <- matrix(vector,byrow = TRUE,nrow = 2)
V
## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 4 3 2 1
dim(vector)
## NULL
dim(vector)<- c(4,2)
vector
## [,1] [,2]
## [1,] 1 4
## [2,] 2 3
## [3,] 3 2
## [4,] 4 1
vector <- t(vector)#matrix transpuesta
vector
## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 4 3 2 1
# Una matrix 4 × 5 de enteros aleatoriamente
# generados por Poisson con promedio 1.5:
X <- matrix(rpois(20,1.5),nrow = 4)
X
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 0 0 3
## [2,] 1 3 2 0 3
## [3,] 0 4 0 1 0
## [4,] 1 2 0 2 1
rownames(X) <- rownames(X,do.NULL = FALSE,prefix = "Trial.")
X
## [,1] [,2] [,3] [,4] [,5]
## Trial.1 0 1 0 0 3
## Trial.2 1 3 2 0 3
## Trial.3 0 4 0 1 0
## Trial.4 1 2 0 2 1
drugname <- c("aspirin","paracetamol","nurofen",
"hedex","placebo")
colnames(X) <- drugname
X
## aspirin paracetamol nurofen hedex placebo
## Trial.1 0 1 0 0 3
## Trial.2 1 3 2 0 3
## Trial.3 0 4 0 1 0
## Trial.4 1 2 0 2 1
mean(X[,5])#Aqui esta el promedio de la columna 5,
## [1] 1.75
#calculada sobre todas las files(en blanco)
var(X[4,])
## [1] 0.7
# Hay funciones especiales para calcular estadisticas
# en matrices:
rowSums(X)
## Trial.1 Trial.2 Trial.3 Trial.4
## 4 9 5 6
colSums(X)
## aspirin paracetamol nurofen hedex placebo
## 2 10 2 3 7
rowMeans(X)
## Trial.1 Trial.2 Trial.3 Trial.4
## 0.8 1.8 1.0 1.2
colMeans(X)
## aspirin paracetamol nurofen hedex placebo
## 0.50 2.50 0.50 0.75 1.75
apply(X,2,mean)
## aspirin paracetamol nurofen hedex placebo
## 0.50 2.50 0.50 0.75 1.75
# para sumar grupos entre las columnas usamos
# rowsum (singular y minuscula)
X
## aspirin paracetamol nurofen hedex placebo
## Trial.1 0 1 0 0 3
## Trial.2 1 3 2 0 3
## Trial.3 0 4 0 1 0
## Trial.4 1 2 0 2 1
# queremos agrupar
# fila 1 y fila 4 (como grupo A) y
# fila 2 y fila 3 (grupo B)
# el vector de agrupacion vector debe tener
# una longitud igual al numero de filas:
grupo <- c("A","B","B","A")
rowsum(X,grupo)
## aspirin paracetamol nurofen hedex placebo
## A 1 3 0 2 4
## B 1 7 2 1 3
#2.8.6 aplicando funciones con apply, sapply and lapply
(X <- matrix(1:24,nrow = 4))
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 5 9 13 17 21
## [2,] 2 6 10 14 18 22
## [3,] 3 7 11 15 19 23
## [4,] 4 8 12 16 20 24
apply(X,1,sum)
## [1] 66 72 78 84
apply(X,2,sum)
## [1] 10 26 42 58 74 90
apply(X,1,sqrt)
## [,1] [,2] [,3] [,4]
## [1,] 1.000000 1.414214 1.732051 2.000000
## [2,] 2.236068 2.449490 2.645751 2.828427
## [3,] 3.000000 3.162278 3.316625 3.464102
## [4,] 3.605551 3.741657 3.872983 4.000000
## [5,] 4.123106 4.242641 4.358899 4.472136
## [6,] 4.582576 4.690416 4.795832 4.898979
apply(X,2,sqrt)
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.000000 2.236068 3.000000 3.605551 4.123106 4.582576
## [2,] 1.414214 2.449490 3.162278 3.741657 4.242641 4.690416
## [3,] 1.732051 2.645751 3.316625 3.872983 4.358899 4.795832
## [4,] 2.000000 2.828427 3.464102 4.000000 4.472136 4.898979
# la matrix resultante es la transpuesta
# Puede proporcionar su propia function definida
# (aqui x2 + x) y aplicarla asi:
apply(X, 1,function(X)X^2+X)#la X es mayuscula
## [,1] [,2] [,3] [,4]
## [1,] 2 6 12 20
## [2,] 30 42 56 72
## [3,] 90 110 132 156
## [4,] 182 210 240 272
## [5,] 306 342 380 420
## [6,] 462 506 552 600
# sapply(3:7,seq)La funcion sapply es muy util
# con complicados calculos iterativos.
# La siguiente data muestra la reaccion en cadena
# en un periodo de 50 dias, usando y = exp(–ax):
sapdecay <- read.table("C:\\Users\\Luis\\Documents\\therbook\\sapdecay.txt",
header = TRUE)
names(sapdecay)
## [1] "x" "y"
head(sapdecay)
## x y
## 1 0 1.0000000
## 2 2 0.9602354
## 3 4 0.8446638
## 4 6 0.7069363
## 5 8 0.7086414
## 6 10 0.6097954
sapdecay
## x y
## 1 0 1.00000000
## 2 2 0.96023540
## 3 4 0.84466377
## 4 6 0.70693633
## 5 8 0.70864140
## 6 10 0.60979536
## 7 12 0.51325264
## 8 14 0.47139969
## 9 16 0.38408127
## 10 18 0.37227780
## 11 20 0.33408743
## 12 22 0.27330733
## 13 24 0.25162333
## 14 26 0.22278651
## 15 28 0.21021653
## 16 30 0.19059182
## 17 32 0.16690720
## 18 34 0.14073164
## 19 36 0.12028785
## 20 38 0.11646071
## 21 40 0.09273687
## 22 42 0.08397261
## 23 44 0.08829140
## 24 46 0.07179059
## 25 48 0.05932741
## 26 50 0.06170623
# escribimos una funcion para calcular
# la suma de los cuadrados de las diferencias
# entre el valor observado de (y)
# y los valores predicciones (yf) de y, when
# dado un valor del parametro amount:
sumsq <- function(a,xv=sapdecay$x,yv=sapdecay$y)
{yf <- exp(-a*xv)
sum((yv-yf)^2)}
# Podemos darnos una idea de la constante de decay, a,
# para esta data mediante la regresion linear log(y)
# vs x, asi
modelD <- lm(log(sapdecay$y)~sapdecay$x)
summary(modelD)
##
## Call:
## lm(formula = log(sapdecay$y) ~ sapdecay$x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08542 -0.04586 0.01362 0.03071 0.09940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0468837 0.0201983 2.321 0.0291 *
## sapdecay$x -0.0584863 0.0006928 -84.421 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05299 on 24 degrees of freedom
## Multiple R-squared: 0.9966, Adjusted R-squared: 0.9965
## F-statistic: 7127 on 1 and 24 DF, p-value: < 2.2e-16
# nuestro parametro es cercano a
# 0.058.Generamos un rango de valores
# con un intervalo a cada lado de 0.058:
a <- seq(0.01,0.2,0.005)
a
## [1] 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060
## [12] 0.065 0.070 0.075 0.080 0.085 0.090 0.095 0.100 0.105 0.110 0.115
## [23] 0.120 0.125 0.130 0.135 0.140 0.145 0.150 0.155 0.160 0.165 0.170
## [34] 0.175 0.180 0.185 0.190 0.195 0.200
plot(a,sapply(a,sumsq),type = "l",
main = "Modelo de desintegracion radioactiva")
a[min(sapply(a,sumsq))==sapply(a,sumsq)]
## [1] 0.055