set.seed(123)
rto = sort.int(rnorm(30, 3, 0.5), partial = 15)
R = list()
r = list()
prom.arit = list()
prom.geo = list()
prom.arm = list()
prom.heroninano = list()
k = list()
m.cov = list()
for (j in 1:15) {
set.seed(123)
pH = sort.int(rnorm(30, 5.5, 0.5), partial = j)
A = sort.int(rnorm(30, 40, 0.5), partial = j)
A = A[length(A):1]
MO = sort.int(rnorm(30, 3, 0.5), partial = j)
M = matrix(cbind(pH,A,MO),ncol = 3)
r[[j]] = abs(cor(M))
R[[j]] = abs(cor(M)[upper.tri(cor(M))])
prom.arit[[j]] = sum(R[[j]])/length(R[[j]])
prom.geo[[j]] = (prod(R[[j]]))^(1/length(R[[j]]))
prom.arm[[j]] = length(R[[j]])/sum(1/R[[j]])
prom.heroninano[[j]] = (prom.arit[[j]]+prom.geo[[j]])/2
# M.e = scale(M)
M.e = M
s.m.e = cov(M.e)
k[[j]] = max(eigen(s.m.e)$values)/min(eigen(s.m.e)$values)
}
r
## [[1]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.4353285 0.6052423
## [2,] 0.4353285 1.0000000 0.5830460
## [3,] 0.6052423 0.5830460 1.0000000
##
## [[2]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.1132112 0.4440392
## [2,] 0.1132112 1.0000000 0.2257159
## [3,] 0.4440392 0.2257159 1.0000000
##
## [[3]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8846405 0.5116056
## [2,] 0.8846405 1.0000000 0.5077232
## [3,] 0.5116056 0.5077232 1.0000000
##
## [[4]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.7380885 0.4846949
## [2,] 0.7380885 1.0000000 0.4754659
## [3,] 0.4846949 0.4754659 1.0000000
##
## [[5]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.7762204 0.3938064
## [2,] 0.7762204 1.0000000 0.4750257
## [3,] 0.3938064 0.4750257 1.0000000
##
## [[6]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8475606 0.7266952
## [2,] 0.8475606 1.0000000 0.7090885
## [3,] 0.7266952 0.7090885 1.0000000
##
## [[7]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8187133 0.8031616
## [2,] 0.8187133 1.0000000 0.6804810
## [3,] 0.8031616 0.6804810 1.0000000
##
## [[8]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.5609480 0.7536904
## [2,] 0.5609480 1.0000000 0.3158857
## [3,] 0.7536904 0.3158857 1.0000000
##
## [[9]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.5096722 0.8014837
## [2,] 0.5096722 1.0000000 0.5060652
## [3,] 0.8014837 0.5060652 1.0000000
##
## [[10]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.6891555 0.6960207
## [2,] 0.6891555 1.0000000 0.7176583
## [3,] 0.6960207 0.7176583 1.0000000
##
## [[11]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8574291 0.6742008
## [2,] 0.8574291 1.0000000 0.7291865
## [3,] 0.6742008 0.7291865 1.0000000
##
## [[12]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.5792782 0.8294224
## [2,] 0.5792782 1.0000000 0.5805707
## [3,] 0.8294224 0.5805707 1.0000000
##
## [[13]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.8669882 0.8486837
## [2,] 0.8669882 1.0000000 0.8505854
## [3,] 0.8486837 0.8505854 1.0000000
##
## [[14]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.6618105 0.6942600
## [2,] 0.6618105 1.0000000 0.8355502
## [3,] 0.6942600 0.8355502 1.0000000
##
## [[15]]
## [,1] [,2] [,3]
## [1,] 1.0000000 0.7672636 0.8236728
## [2,] 0.7672636 1.0000000 0.8433691
## [3,] 0.8236728 0.8433691 1.0000000
R
## [[1]]
## [1] 0.4353285 0.6052423 0.5830460
##
## [[2]]
## [1] 0.1132112 0.4440392 0.2257159
##
## [[3]]
## [1] 0.8846405 0.5116056 0.5077232
##
## [[4]]
## [1] 0.7380885 0.4846949 0.4754659
##
## [[5]]
## [1] 0.7762204 0.3938064 0.4750257
##
## [[6]]
## [1] 0.8475606 0.7266952 0.7090885
##
## [[7]]
## [1] 0.8187133 0.8031616 0.6804810
##
## [[8]]
## [1] 0.5609480 0.7536904 0.3158857
##
## [[9]]
## [1] 0.5096722 0.8014837 0.5060652
##
## [[10]]
## [1] 0.6891555 0.6960207 0.7176583
##
## [[11]]
## [1] 0.8574291 0.6742008 0.7291865
##
## [[12]]
## [1] 0.5792782 0.8294224 0.5805707
##
## [[13]]
## [1] 0.8669882 0.8486837 0.8505854
##
## [[14]]
## [1] 0.6618105 0.6942600 0.8355502
##
## [[15]]
## [1] 0.7672636 0.8236728 0.8433691
prom.arit[[1]]
## [1] 0.5412056
prom.geo[[1]]
## [1] 0.5355702
prom.arm[[1]]
## [1] 0.5296164
prom.heroninano[[1]]
## [1] 0.5383879
plot(unlist(prom.arit),unlist(k))
########
# Graficar diagrama de venn para las varaibles
# explicativas NPK representar la correlacion
# Libro gujarati
set.seed(123)
pH = sort.int(rnorm(30, 5.5, 0.5), partial = 15)
rto = sort.int(rnorm(30, 3.0, 0.5), partial = 15)
pH2 = rep(c(4,4.5,5,5.5,6), each = 6)
plot(pH,rto, pch = 16)
grid(10, col = gray(0.5))
mod1 = lm(rto~pH)
abline(mod1)
v.pronos = mod1$fitted.values
points(pH, v.pronos, col ='green', pch = 16)
plot(pH2, rto, pch = 16)
summary(mod1)
##
## Call:
## lm(formula = rto ~ pH)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68428 -0.06695 -0.02683 0.06663 0.55015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.73545 0.50573 -1.454 0.157
## pH 0.69838 0.09199 7.592 2.86e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.243 on 28 degrees of freedom
## Multiple R-squared: 0.673, Adjusted R-squared: 0.6614
## F-statistic: 57.64 on 1 and 28 DF, p-value: 2.862e-08
# Ajustando modelo por minimos cuadrados
(mod1 = lm(rto~pH))
##
## Call:
## lm(formula = rto ~ pH)
##
## Coefficients:
## (Intercept) pH
## -0.7354 0.6984
(mod2 = lm(rto~pH2))
##
## Call:
## lm(formula = rto ~ pH2)
##
## Coefficients:
## (Intercept) pH2
## 0.4775 0.5223
e1 = mod1$residuals
plot(e1, type ='l')
e2 = mod2$residuals
plot(e2)
hist(e1)
hist(e2)
# Supuesto 2
cor(e1,pH)
## [1] -9.062801e-16
cor(e2,pH2)
## [1] -4.546989e-16
# Supuesto 3
mean(e1)
## [1] 1.01644e-17
mean(e2)
## [1] 4.681043e-18
mean(e2[pH2==4.0])
## [1] -0.01100354
mean(e2[pH2==4.5])
## [1] 0.04646977
mean(e2[pH2==5.0])
## [1] -0.03945873
mean(e2[pH2==5.5])
## [1] -0.01647767
mean(e2[pH2==6.0])
## [1] 0.02047018
# Supuesto 4
var(e1)
## [1] 0.05700962
var(e2)
## [1] 0.03324372
boxplot(e2~pH2)
tapply(e2, pH2, var)
## 4 4.5 5 5.5 6
## 0.062928042 0.006600074 0.018155683 0.012980947 0.086715151
bartlett.test(e2~pH2)
##
## Bartlett test of homogeneity of variances
##
## data: e2 by pH2
## Bartlett's K-squared = 10.02, df = 4, p-value = 0.04009
# Supuesto 5
# Los residuales no muestran autocorrelacion
acf(e1)
acf(e2)
set.seed(123)
(mes = rep(1:12, 2))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11
## [24] 12
(pcp = runif(24, 40, 60))
## [1] 45.75155 55.76610 48.17954 57.66035 58.80935 40.91113 50.56211
## [8] 57.84838 51.02870 49.13229 59.13667 49.06668 53.55141 51.45267
## [15] 42.05849 57.99650 44.92175 40.84119 46.55841 59.09007 57.79079
## [22] 53.85607 52.81014 59.88540
pcp[mes == 5] = pcp[mes == 5]*2
pcp[mes == 10] = pcp[mes == 10]*2.2
acf(pcp)
# Supuesto 7
cic = 2*pH
plot(cic,pH)
(mod3 = lm(rto ~ pH + cic))
##
## Call:
## lm(formula = rto ~ pH + cic)
##
## Coefficients:
## (Intercept) pH cic
## -0.7354 0.6984 NA
set.seed(123)
A = round(rnorm(30, 40), 1)
Ar = round(rnorm(30, 20), 1)
L = 100 - (A+Ar)
(mod4 = lm(rto ~ pH + A+Ar+L))
##
## Call:
## lm(formula = rto ~ pH + A + Ar + L)
##
## Coefficients:
## (Intercept) pH A Ar L
## -0.20520 0.70073 0.01430 -0.05522 NA
cor(A,L)
## [1] -0.7220848
cor(A,Ar)
## [1] -0.146547
# Supuesto 10
shapiro.test(e1)
##
## Shapiro-Wilk normality test
##
## data: e1
## W = 0.9126, p-value = 0.01732
shapiro.test(e2)
##
## Shapiro-Wilk normality test
##
## data: e2
## W = 0.98557, p-value = 0.9465
#####
# summary(aov(rto ~ pH))
# pH, A
MO = rnorm(30, 3, 0.6)
(X = matrix(c(pH,A,MO), ncol = 3))
## [,1] [,2] [,3]
## [1,] 5.219762 39.4 3.227784
## [2,] 4.930932 39.8 2.698606
## [3,] 4.656653 41.6 2.800076
## [4,] 5.187480 40.1 2.388855
## [5,] 5.135554 40.1 2.356925
## [6,] 4.986998 41.7 3.182117
## [7,] 4.966088 40.5 3.268926
## [8,] 4.867469 38.7 3.031803
## [9,] 5.156574 39.3 3.553360
## [10,] 4.516691 39.6 4.230051
## [11,] 5.222079 41.2 2.705381
## [12,] 5.263604 40.4 1.614499
## [13,] 5.277169 40.4 3.603443
## [14,] 5.384911 40.1 2.574480
## [15,] 5.391013 39.4 2.587195
## [16,] 5.535254 41.8 3.615343
## [17,] 5.748925 40.5 2.829136
## [18,] 5.700386 38.0 2.267569
## [19,] 5.850678 40.7 3.108782
## [20,] 5.679907 39.5 2.916665
## [21,] 5.730458 38.9 3.003459
## [22,] 5.576687 39.8 3.231168
## [23,] 5.918894 39.0 2.777604
## [24,] 5.564644 39.3 3.386626
## [25,] 5.555341 39.4 2.867708
## [26,] 6.279354 38.3 3.199069
## [27,] 6.357532 40.8 3.658103
## [28,] 6.393457 40.2 3.261109
## [29,] 6.112041 38.9 2.804441
## [30,] 6.126907 41.3 3.689285
cor(X)
## [,1] [,2] [,3]
## [1,] 1.00000000 -0.1745336 0.06573727
## [2,] -0.17453361 1.0000000 0.16690999
## [3,] 0.06573727 0.1669100 1.00000000
s = cov(X)
val.p = eigen(s)$values
vec.p = eigen(s)$vectors
val.p; vec.p
## [1] 0.9963006 0.2775038 0.2161787
## [,1] [,2] [,3]
## [1,] 0.1080182 0.49988416 0.8593299
## [2,] -0.9874830 -0.04598009 0.1508744
## [3,] -0.1149318 0.86487089 -0.4886604
diag(val.p)
## [,1] [,2] [,3]
## [1,] 0.9963006 0.0000000 0.0000000
## [2,] 0.0000000 0.2775038 0.0000000
## [3,] 0.0000000 0.0000000 0.2161787
vec.p %*% diag(val.p) %*% t(vec.p)
## [,1] [,2] [,3]
## [1,] 0.24060531 -0.08462215 0.01682797
## [2,] -0.08462215 0.97702299 0.08609983
## [3,] 0.01682797 0.08609983 0.27235485
s
## [,1] [,2] [,3]
## [1,] 0.24060531 -0.08462215 0.01682797
## [2,] -0.08462215 0.97702299 0.08609983
## [3,] 0.01682797 0.08609983 0.27235485
(k = max(val.p)/min(val.p))
## [1] 4.608689
# Ordenados
(X = matrix(c(sort(pH),sort(A),sort(MO)), ncol = 3))
## [,1] [,2] [,3]
## [1,] 4.516691 38.0 1.614499
## [2,] 4.656653 38.3 2.267569
## [3,] 4.867469 38.7 2.356925
## [4,] 4.930932 38.9 2.388855
## [5,] 4.966088 38.9 2.574480
## [6,] 4.986998 39.0 2.587195
## [7,] 5.135554 39.3 2.698606
## [8,] 5.156574 39.3 2.705381
## [9,] 5.187480 39.4 2.777604
## [10,] 5.219762 39.4 2.800076
## [11,] 5.222079 39.4 2.804441
## [12,] 5.263604 39.5 2.829136
## [13,] 5.277169 39.6 2.867708
## [14,] 5.384911 39.8 2.916665
## [15,] 5.391013 39.8 3.003459
## [16,] 5.535254 40.1 3.031803
## [17,] 5.555341 40.1 3.108782
## [18,] 5.564644 40.1 3.182117
## [19,] 5.576687 40.2 3.199069
## [20,] 5.679907 40.4 3.227784
## [21,] 5.700386 40.4 3.231168
## [22,] 5.730458 40.5 3.261109
## [23,] 5.748925 40.5 3.268926
## [24,] 5.850678 40.7 3.386626
## [25,] 5.918894 40.8 3.553360
## [26,] 6.112041 41.2 3.603443
## [27,] 6.126907 41.3 3.615343
## [28,] 6.279354 41.6 3.658103
## [29,] 6.357532 41.7 3.689285
## [30,] 6.393457 41.8 4.230051
cor(X)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9995329 0.9736389
## [2,] 0.9995329 1.0000000 0.9735269
## [3,] 0.9736389 0.9735269 1.0000000
(s = cov(X))
## [,1] [,2] [,3]
## [1,] 0.2406053 0.4846208 0.2492402
## [2,] 0.4846208 0.9770230 0.5021898
## [3,] 0.2492402 0.5021898 0.2723548
val.p = eigen(s)$values
vec.p = eigen(s)$vectors
(k = max(val.p)/min(val.p))
## [1] 8221.33
# Estandarizados
pH.e = (pH - mean(pH))/sd(pH)
A.e = (A - mean(A))/sd(A)
MO.e = (MO - mean(MO))/sd(MO)
(X = matrix(c(sort(pH.e),sort(A.e),sort(MO.e)), ncol = 3))
## [,1] [,2] [,3]
## [1,] -1.9566293 -1.9795408 -2.68292242
## [2,] -1.6712928 -1.6760337 -1.43153252
## [3,] -1.2415080 -1.2713575 -1.26031207
## [4,] -1.1121295 -1.0690195 -1.19912993
## [5,] -1.0404567 -1.0690195 -0.84344262
## [6,] -0.9978288 -0.9678504 -0.81907805
## [7,] -0.6949706 -0.6643433 -0.60559623
## [8,] -0.6521193 -0.6643433 -0.59261352
## [9,] -0.5891105 -0.5631743 -0.45422309
## [10,] -0.5232985 -0.5631743 -0.41116386
## [11,] -0.5185744 -0.5631743 -0.40279889
## [12,] -0.4339188 -0.4620053 -0.35547897
## [13,] -0.4062648 -0.3608362 -0.28156898
## [14,] -0.1866137 -0.1584982 -0.18775916
## [15,] -0.1741751 -0.1584982 -0.02144900
## [16,] 0.1198863 0.1450090 0.03286277
## [17,] 0.1608374 0.1450090 0.18036814
## [18,] 0.1798022 0.1450090 0.32089014
## [19,] 0.2043533 0.2461780 0.35337292
## [20,] 0.4147858 0.4485161 0.40839460
## [21,] 0.4565354 0.4485161 0.41487995
## [22,] 0.5178431 0.5496851 0.47225113
## [23,] 0.5554915 0.5496851 0.48722972
## [24,] 0.7629319 0.7520232 0.71276222
## [25,] 0.9020011 0.8531922 1.03225280
## [26,] 1.2957653 1.2578683 1.12821929
## [27,] 1.3260734 1.3590374 1.15102106
## [28,] 1.6368622 1.6625445 1.23295732
## [29,] 1.7962422 1.7637135 1.29270551
## [30,] 1.8694796 1.8648826 2.32890176
cor(X)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9995329 0.9736389
## [2,] 0.9995329 1.0000000 0.9735269
## [3,] 0.9736389 0.9735269 1.0000000
(s = cov(X))
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9995329 0.9736389
## [2,] 0.9995329 1.0000000 0.9735269
## [3,] 0.9736389 0.9735269 1.0000000
val.p = eigen(s)$values
vec.p = eigen(s)$vectors
val.p; vec.p
## [1] 2.9645168764 0.0350161113 0.0004670124
## [,1] [,2] [,3]
## [1,] 0.5790608 -0.4045052 0.707858790
## [2,] 0.5790392 -0.4071610 -0.706352279
## [3,] 0.5739357 0.8188989 -0.001546665
diag(val.p)
## [,1] [,2] [,3]
## [1,] 2.964517 0.00000000 0.0000000000
## [2,] 0.000000 0.03501611 0.0000000000
## [3,] 0.000000 0.00000000 0.0004670124
vec.p %*% diag(val.p) %*% t(vec.p)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9995329 0.9736389
## [2,] 0.9995329 1.0000000 0.9735269
## [3,] 0.9736389 0.9735269 1.0000000
s
## [,1] [,2] [,3]
## [1,] 1.0000000 0.9995329 0.9736389
## [2,] 0.9995329 1.0000000 0.9735269
## [3,] 0.9736389 0.9735269 1.0000000
(k = max(val.p)/min(val.p))
## [1] 6347.834