CORRELACION CANONICA
Melissa Ruiz
20/6/2021
R Markdown
##This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com. ## ##When you click the Knit button a document will be generated that includes both content as well as ##the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
##{r cars} ##summary(cars) ##
Including Plots
##You can also embed plots, for example:
##{r pressure, echo=FALSE} ##plot(pressure) ##
##Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code ##that generated the plot.
data<-read.table("BD_CANONICA.csv",sep=";")
colnames(data) <-c("Departamento","Humedad [%]","Temperatura promedio [C]","Precipitación [mm]","Producción de Agua [miles de m3]","Caudal Promedio [m3/s]")
data## Departamento Humedad [%] Temperatura promedio [C] Precipitación [mm]
## 1 AHUACHAPÁN 73.08 23.5 136
## 2 SANTA ANA 68.50 24.5 134
## 3 SONSONATE 68.30 26.5 144
## 4 CHALATENANGO 74.20 24.0 146
## 5 LA LIBERTAD 74.80 28.0 161
## 6 SAN SALVADOR 73.30 24.0 146
## 7 CUSCATLÁN 76.60 25.0 142
## 8 LA PAZ 75.00 27.5 147
## 9 CABAÑAS 67.00 25.5 139
## 10 SAN VICENTE 69.80 25.5 142
## 11 USULUTÁN 80.00 28.0 155
## 12 SAN MIGUEL 58.00 29.0 157
## 13 MORAZÁN 76.00 27.0 151
## 14 LA UNIÓN 68.00 30.0 165
## Producción de Agua [miles de m3] Caudal Promedio [m3/s]
## 1 12666.5 25.60
## 2 49344.9 57.03
## 3 17259.4 2.96
## 4 3335.1 1.93
## 5 31390.6 10.41
## 6 177522.5 7.86
## 7 12592.4 1.92
## 8 12088.5 2.94
## 9 5520.0 2.94
## 10 8366.0 7.86
## 11 14627.8 2.94
## 12 15266.2 19.01
## 13 488.3 30.65
## 14 10655.4 12.59cov(data[,-1])## Humedad [%] Temperatura promedio [C]
## Humedad [%] 30.106007 -2.588571
## Temperatura promedio [C] -2.588571 4.104396
## Precipitación [mm] -1.605385 16.076923
## Producción de Agua [miles de m3] 12124.110066 -29635.520330
## Caudal Promedio [m3/s] -19.366204 -4.549286
## Precipitación [mm]
## Humedad [%] -1.605385
## Temperatura promedio [C] 16.076923
## Precipitación [mm] 85.500000
## Producción de Agua [miles de m3] -34235.507692
## Caudal Promedio [m3/s] -38.370769
## Producción de Agua [miles de m3]
## Humedad [%] 12124.11
## Temperatura promedio [C] -29635.52
## Precipitación [mm] -34235.51
## Producción de Agua [miles de m3] 2040914703.80
## Caudal Promedio [m3/s] 45727.62
## Caudal Promedio [m3/s]
## Humedad [%] -19.366204
## Temperatura promedio [C] -4.549286
## Precipitación [mm] -38.370769
## Producción de Agua [miles de m3] 45727.621418
## Caudal Promedio [m3/s] 242.534213##la matriz de covarianzas de X, Sxx es:
X <-data[,(2:4)]
cov(X)## Humedad [%] Temperatura promedio [C]
## Humedad [%] 30.106007 -2.588571
## Temperatura promedio [C] -2.588571 4.104396
## Precipitación [mm] -1.605385 16.076923
## Precipitación [mm]
## Humedad [%] -1.605385
## Temperatura promedio [C] 16.076923
## Precipitación [mm] 85.500000##La matriz de covarianzas de Y,Syy es:
Y <-data[,(5:6)]
cov(Y)## Producción de Agua [miles de m3]
## Producción de Agua [miles de m3] 2.040915e+09
## Caudal Promedio [m3/s] 4.572762e+04
## Caudal Promedio [m3/s]
## Producción de Agua [miles de m3] 45727.6214
## Caudal Promedio [m3/s] 242.5342##La matriz de covarianzas de X e Y, Sxy es:
cov(X,Y)## Producción de Agua [miles de m3]
## Humedad [%] 12124.11
## Temperatura promedio [C] -29635.52
## Precipitación [mm] -34235.51
## Caudal Promedio [m3/s]
## Humedad [%] -19.366204
## Temperatura promedio [C] -4.549286
## Precipitación [mm] -38.370769##encontrar los autovalores y autovectores de las siguientes matrices cuadradas
A <-solve(cov(X,X))%*%cov(X,Y)%*%solve(cov(Y,Y))%*%cov(Y,X)
B <-solve(cov(Y,Y))%*%cov(Y,X)%*%solve(cov(X,X))%*%cov(X,Y)
A## Humedad [%] Temperatura promedio [C]
## Humedad [%] 0.03499497 0.03028919
## Temperatura promedio [C] -0.29296465 0.32816096
## Precipitación [mm] 0.08912795 -0.04824701
## Precipitación [mm]
## Humedad [%] 0.10891510
## Temperatura promedio [C] 0.13150253
## Precipitación [mm] 0.05251439B## Producción de Agua [miles de m3]
## Producción de Agua [miles de m3] 0.2810228
## Caudal Promedio [m3/s] -136.7304932
## Caudal Promedio [m3/s]
## Producción de Agua [miles de m3] -1.296891e-05
## Caudal Promedio [m3/s] 1.346475e-01##Los autovalores son las soluciones de la ecuación |A–λI|=0 . En este caso tenemos tres autovalores ( el mínimo entre p=3 y q=2 es 2 ):
eigen(A)$values## [1] 2.922726e-01 1.233978e-01 -6.938894e-18eigen(B)$values ## [1] 0.2922726 0.1233978##Las correspondientes correlaciones canónicas son las raíces cuadradas positivas de los autovalores anteriores:
r <-sqrt(abs(eigen(A)$values))
r## [1] 5.406224e-01 3.512802e-01 2.634178e-09##Los primeros vectores canónicos se obtienen a partir del primer autovector de A y del primer autovector de B:
a <-eigen(A)$vector[,1]
b <-eigen(B)$vector[,1]
a## [1] 0.03793359 0.98227561 -0.18356377b## [1] 0.001152815 -0.999999336##Como atSxxa debe ser igual a 1 y
t(a)%*%cov(X,X)%*%a## [,1]
## [1,] 0.9162729t(b)%*%cov(Y,Y)%*%b## [,1]
## [1,] 2849.442##vectores canonicos
a1<-(a*(1/(1.54201^0.5)))
a1## [1] 0.03054783 0.79102420 -0.14782347b1<-(b*(1/(1.065463^0.5)))
b1## [1] 0.001116838 -0.968791958##Estimando los segundos vectores canonicos
a11 <-eigen(A)$vector[,2]
b11 <-eigen(B)$vector[,2]
a11## [1] 0.6350356 0.7031329 0.3198968b11## [1] 8.227694e-05 1.000000e+00##Como atSxxa debe ser igual a 1 y
t(a11)%*%cov(X,X)%*%a11## [,1]
## [1,] 27.18804t(b11)%*%cov(Y,Y)%*%b11## [,1]
## [1,] 263.8748##vectores canonicos
a2<-(a*(1/(1.54201^0.5)))
a2## [1] 0.03054783 0.79102420 -0.14782347b2<-(b*(1/(1.065463^0.5)))
b2## [1] 0.001116838 -0.968791958#TEST PARA DESARROLLO DE EJERCICIO CANONICA
library(GGally)## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(fields)## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.6-0 (2020-12-14) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: viridis## Loading required package: viridisLite## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source codelibrary(CCA)## Loading required package: fda## Loading required package: splines## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:spam':
##
## det## Loading required package: fds## Loading required package: rainbow## Loading required package: MASS## Loading required package: pcaPP## Loading required package: RCurl##
## Attaching package: 'fda'## The following object is masked from 'package:graphics':
##
## matplotlibrary(vegan)## Loading required package: permute## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:fda':
##
## melanoma## This is vegan 2.5-7clima<-read.table("CLIMA.csv",sep=";")
agua<-read.table("AGUA.csv",sep=";")
colnames(agua)<-c("Produccion de agua [miles de m3]","Caudal Promedio [m3/s]")
colnames(clima)<-c("Humedad [%]","Temperatura [C]","Precipitacion [mm]")
clima## Humedad [%] Temperatura [C] Precipitacion [mm]
## 1 73.08 23.5 136
## 2 68.50 24.5 134
## 3 68.30 26.5 144
## 4 74.20 24.0 146
## 5 74.80 28.0 161
## 6 73.30 24.0 146
## 7 76.60 25.0 142
## 8 75.00 27.5 147
## 9 67.00 25.5 139
## 10 69.80 25.5 142
## 11 80.00 28.0 155
## 12 58.00 29.0 157
## 13 76.00 27.0 151
## 14 68.00 30.0 165agua## Produccion de agua [miles de m3] Caudal Promedio [m3/s]
## 1 12666.5 25.60
## 2 49344.9 57.03
## 3 17259.4 2.96
## 4 3335.1 1.93
## 5 31390.6 10.41
## 6 177522.5 7.86
## 7 12592.4 1.92
## 8 12088.5 2.94
## 9 5520.0 2.94
## 10 8366.0 7.86
## 11 14627.8 2.94
## 12 15266.2 19.01
## 13 488.3 30.65
## 14 10655.4 12.59#gráficas la correlación dentro de cada conjunto de variables
ggpairs(agua, title = "Disponibilidad de Agua")
ggpairs(clima, title = "Variables climaticas")
#Para ambas variables
deptos <-cbind(agua,clima)
ggduo(deptos,columnsX = 1:2,columnsY = 3:5,
types = list(continuous = "smooth_lm"),
title = "Correlación entre variables Disponibilidad de agua y climaticas",
xlab = "Disponibilidad de agua",
ylab = "Clima"
)
#Correlacion canónica
mat_cor <-matcor(clima,agua)
mat_cor## $Xcor
## Humedad [%] Temperatura [C] Precipitacion [mm]
## Humedad [%] 1.0000000 -0.2328675 -0.0316424
## Temperatura [C] -0.2328675 1.0000000 0.8582134
## Precipitacion [mm] -0.0316424 0.8582134 1.0000000
##
## $Ycor
## Produccion de agua [miles de m3]
## Produccion de agua [miles de m3] 1.00000000
## Caudal Promedio [m3/s] 0.06499496
## Caudal Promedio [m3/s]
## Produccion de agua [miles de m3] 0.06499496
## Caudal Promedio [m3/s] 1.00000000
##
## $XYcor
## Humedad [%] Temperatura [C] Precipitacion [mm]
## Humedad [%] 1.00000000 -0.2328675 -0.0316424
## Temperatura [C] -0.23286755 1.0000000 0.8582134
## Precipitacion [mm] -0.03164240 0.8582134 1.0000000
## Produccion de agua [miles de m3] 0.04891149 -0.3237990 -0.0819561
## Caudal Promedio [m3/s] -0.22663732 -0.1441890 -0.2664593
## Produccion de agua [miles de m3]
## Humedad [%] 0.04891149
## Temperatura [C] -0.32379901
## Precipitacion [mm] -0.08195610
## Produccion de agua [miles de m3] 1.00000000
## Caudal Promedio [m3/s] 0.06499496
## Caudal Promedio [m3/s]
## Humedad [%] -0.22663732
## Temperatura [C] -0.14418899
## Precipitacion [mm] -0.26645934
## Produccion de agua [miles de m3] 0.06499496
## Caudal Promedio [m3/s] 1.00000000cca1<- cc(clima,agua)
cca1## $cor
## [1] 0.5406224 0.3512802
##
## $names
## $names$Xnames
## [1] "Humedad [%]" "Temperatura [C]" "Precipitacion [mm]"
##
## $names$Ynames
## [1] "Produccion de agua [miles de m3]" "Caudal Promedio [m3/s]"
##
## $names$ind.names
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"
##
##
## $xcoef
## [,1] [,2]
## Humedad [%] 0.03962886 -0.12178929
## Temperatura [C] 1.02617386 -0.13484922
## Precipitacion [mm] -0.19176731 -0.06135089
##
## $ycoef
## [,1] [,2]
## Produccion de agua [miles de m3] -2.159631e-05 5.064997e-06
## Caudal Promedio [m3/s] 1.873354e-02 6.156034e-02
##
## $scores
## $scores$xscores
## [,1] [,2]
## [1,] -0.59516198 0.90250438
## [2,] 0.63304631 1.44815188
## [3,] 0.75979521 0.58930237
## [4,] -1.95536378 0.08516684
## [5,] -0.70340058 -1.44756700
## [6,] -1.99102975 0.19477720
## [7,] -0.06701143 -0.09657309
## [8,] 1.47618053 -0.54558774
## [9,] 0.64094035 1.18923212
## [10,] 0.17659925 0.66416945
## [11,] 0.65327333 -1.71276592
## [12,] 0.42407764 0.70904735
## [13,] 0.23565324 -0.84535599
## [14,] 0.31240167 -1.13450187
##
## $scores$yscores
## [,1] [,2]
## [1,] 0.528777048 0.6851461
## [2,] 0.325453929 2.8057636
## [3,] 0.005460003 -0.6853170
## [4,] 0.286878019 -0.8192507
## [5,] -0.160156968 -0.1551180
## [6,] -3.363837967 0.4280607
## [7,] 0.086767122 -0.7729781
## [8,] 0.116757715 -0.7127388
## [9,] 0.258613107 -0.7460082
## [10,] 0.289319009 -0.4287164
## [11,] 0.061918193 -0.6998773
## [12,] 0.349179084 0.2926309
## [13,] 0.886385659 0.9343433
## [14,] 0.328486048 -0.1259402
##
## $scores$corr.X.xscores
## [,1] [,2]
## Humedad [%] -0.210574015 -0.5866762
## Temperatura [C] 0.506539092 -0.6044372
## Precipitacion [mm] 0.004108974 -0.7806033
##
## $scores$corr.Y.xscores
## [,1] [,2]
## Produccion de agua [miles de m3] -0.5172045 0.1022682
## Caudal Promedio [m3/s] 0.1234430 0.3420003
##
## $scores$corr.X.yscores
## [,1] [,2]
## Humedad [%] -0.113841025 -0.2060877
## Temperatura [C] 0.273846370 -0.2123268
## Precipitacion [mm] 0.002221403 -0.2742105
##
## $scores$corr.Y.yscores
## [,1] [,2]
## Produccion de agua [miles de m3] -0.9566835 0.2911301
## Caudal Promedio [m3/s] 0.2283349 0.9735826#Gráficos de correspondencia canónica
img.matcor(mat_cor, type = 2)
plt.cc(cca1,var.label=T)
cca2<- cca(clima,agua)
cca2## Call: cca(X = clima, Y = agua)
##
## Inertia Proportion Rank
## Total 2.146e-03 1.000e+00
## Constrained 6.573e-05 3.063e-02 2
## Unconstrained 2.080e-03 9.694e-01 2
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2
## 6.537e-05 3.500e-07
##
## Eigenvalues for unconstrained axes:
## CA1 CA2
## 0.0020012 0.0000789plot(cca2)