library(MASS)

sp <-list()
v1<-list()
v2<-list()
v3<-vector()
sd = 2
n = 20
vmG <- 20
vmp <- 5
res <- list()

simular nidos no colineares

F1G <- data.frame(sp = rep("F1G", n))


F1G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F1G$v2 <- rnorm(n = 20, mean = vmG * 2, sd = sd)

F1G$v3 <- rnorm(n = 20, mean = vmG * 4, sd = sd)

res[[length(res)+1]] <- F1G

# triang peq

vm <- 5

F1p <- data.frame(sp = rep("F1p", n))


F1p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F1p$v2 <- rnorm(n = 20, mean = vmp * 2, sd = sd)

F1p$v3 <- rnorm(n = 20, mean = vmp * 4, sd = sd)

res[[length(res)+1]] <- F1p

#F2

#grande

F2G <- data.frame(sp = rep("F2G", n))


F2G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F2G$v2 <- rnorm(n = 20, mean = vmG * 5, sd = sd)

F2G$v3 <- rnorm(n = 20, mean = vmG * 13, sd = sd)

res[[length(res)+1]] <- F2G

# peq


F2p <- data.frame(sp = rep("F2p", n))


F2p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F2p$v2 <- rnorm(n = 20, mean = vmp * 5, sd = sd)

F2p$v3 <- rnorm(n = 20, mean = vmp * 13, sd = sd)

res[[length(res)+1]] <- F2p


#F3

#grande

F3G <- data.frame(sp = rep("F3G", n))


F3G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F3G$v2 <- rnorm(n = 20, mean = vmG * 9, sd = sd)

F3G$v3 <- rnorm(n = 20, mean = vmG * 1.3, sd = sd)

res[[length(res)+1]] <- F3G

# peq


F3p <- data.frame(sp = rep("F3p", n))


F3p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F3p$v2 <- rnorm(n = 20, mean = vmp * 9, sd = sd)

F3p$v3 <- rnorm(n = 20, mean = vmp * 1.3, sd = sd)

res[[length(res)+1]] <- F3p


resT <- do.call(rbind, res)

names(resT)
## [1] "sp" "v1" "v2" "v3"

Analisis discriminante

Datos originales

#### Disc datos originales

ad <- lda(sp ~ v1 + v2 + v3, data = resT, CV = T)

tab1 <- table(resT$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G  20   0   0   0   0   0
##   F1p   0  20   0   0   0   0
##   F2G   0   0  20   0   0   0
##   F2p   0   0   0  20   0   0
##   F3G   0   0   0   0  20   0
##   F3p   0   0   0   0   0  20
diag(prop.table(tab1,1))
## F1G F1p F2G F2p F3G F3p 
##   1   1   1   1   1   1
sum(diag(prop.table(tab1)))
## [1] 1
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(resT$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
## F1 F2 F3 
##  1  1  1
sum(diag(prop.table(tab2)))
## [1] 1

proporcionales al promedio de cada variable dentro de especies

#### Disc proporcionales al promedio de cada variable dentro de especies
dpvs <- lapply(unique(resT$sp), function(x){
y <- resT[resT$sp == x, ] 

y[,-1] <- y[,-1]/colMeans(y[,-1])
return(y)
})

dpvs <- do.call(rbind, dpvs)


ad <- lda(sp ~ v1 + v2 + v3, data = dpvs, CV = T)

tab1 <- table(dpvs$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G   8  12   0   0   0   0
##   F1p  13   4   0   2   1   0
##   F2G   7   0   0   6   0   7
##   F2p   5   5   8   0   0   2
##   F3G   6   0   0   0   4  10
##   F3p   4   2   0   0  13   1
diag(prop.table(tab1,1))
##  F1G  F1p  F2G  F2p  F3G  F3p 
## 0.40 0.20 0.00 0.00 0.20 0.05
sum(diag(prop.table(tab1)))
## [1] 0.1416667
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dpvs$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
##    F1    F2    F3 
## 0.925 0.350 0.700
sum(diag(prop.table(tab2)))
## [1] 0.6583333

proporcionales a una medida dentro de especie

#### Disc proporcionales a una medida dentro de especie
#prop a v1
dv1s <- lapply(unique(resT$sp), function(x){
  y <- resT[resT$sp == x, ] 
  
  y[,-1] <- y[,-1]/y$v2
  return(y)
})

dv1s <- do.call(rbind, dv1s)


ad <- lda(sp ~  v1 + v3, data = dv1s, CV = T)

tab1 <- table(dv1s$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G  16   4   0   0   0   0
##   F1p  12   4   3   1   0   0
##   F2G   0   0  18   2   0   0
##   F2p   0   0   8  12   0   0
##   F3G   0   0   0   0  15   5
##   F3p   0   0   0   0   9  11
diag(prop.table(tab1,1))
##  F1G  F1p  F2G  F2p  F3G  F3p 
## 0.80 0.20 0.90 0.60 0.75 0.55
sum(diag(prop.table(tab1)))
## [1] 0.6333333
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dv1s$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
##  F1  F2  F3 
## 0.9 1.0 1.0
sum(diag(prop.table(tab2)))
## [1] 0.9666667

proporcionales a tamano de cada especie

#### 

dTSs <- lapply(unique(resT$sp), function(x){
  y <- resT[resT$sp == x, ] 
  y[,-1] <- y[,-1]/mean(unlist(y[,-1]))
  return(y)
})

dTSs <- do.call(rbind, dTSs)


ad <- lda(sp ~  v1 + v2 + v3, data = dTSs, CV = T)

# tab1 <- table(dTSs$sp, ad$class)
# 
# tab1
# 
# diag(prop.table(tab1,1))
# 
# sum(diag(prop.table(tab1)))
# 

cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dTSs$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
##    F1    F2    F3 
## 0.925 1.000 1.000
sum(diag(prop.table(tab2)))
## [1] 0.975

NIDOS CON VARIABLES COLINEALES

res2 <- list()
sdcoll <- 0.01

F1G <- data.frame(sp = rep("F1G", n))


F1G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F1G$v2 <- (F1G$v1 * 2) + rnorm(n = 20, mean = 0, sd = sdcoll) * 2

F1G$v3 <- (F1G$v1 * 4) + rnorm(n = 20, mean = 0, sd = sdcoll) * 4 

res2[[length(res2)+1]] <- F1G

# triang peq

vm <- 5

F1p <- data.frame(sp = rep("F1p", n))


F1p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F1p$v2 <- (F1p$v1 * 2) + rnorm(n = 20, mean = 0, sd = sdcoll) * 2

F1p$v3 <- (F1p$v1 * 4) + rnorm(n = 20, mean = 0, sd = sdcoll) * 4

res2[[length(res2)+1]] <- F1p

#F2

#grande

F2G <- data.frame(sp = rep("F2G", n))


F2G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F2G$v2 <- (F2G$v1 * 5) + rnorm(n = 20, mean = 0, sd = sdcoll) * 5

F2G$v3 <- (F2G$v1 * 13) + rnorm(n = 20, mean = 0, sd = sdcoll) * 13

res2[[length(res2)+1]] <- F2G

# peq


F2p <- data.frame(sp = rep("F2p", n))


F2p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F2p$v2 <- (F2p$v1 * 5) + rnorm(n = 20, mean = 0, sd = sdcoll) * 5

F2p$v3 <- (F2p$v1 * 13) + rnorm(n = 20, mean = 0, sd = sdcoll) * 13

res2[[length(res2)+1]] <- F2p


#F3

#grande

F3G <- data.frame(sp = rep("F3G", n))


F3G$v1 <- rnorm(n = 20, mean = vmG, sd = sd)

F3G$v2 <- (F3G$v1 * 9) + rnorm(n = 20, mean = 0, sd = sdcoll) * 9

F3G$v3 <- (F3G$v1 * 1.3) + rnorm(n = 20, mean = 0, sd = sdcoll) *1.3

res2[[length(res2)+1]] <- F3G

# peq


F3p <- data.frame(sp = rep("F3p", n))


F3p$v1 <- rnorm(n = 20, mean = vmp, sd = sd)

F3p$v2 <- (F3p$v1 * 9) + rnorm(n = 20, mean = 0, sd = sdcoll) * 9

F3p$v3 <- (F3p$v1 * 1.3) + rnorm(n = 20, mean = 0, sd = sdcoll) * 1.3


res2[[length(res2)+1]] <- F3p


res2T <- do.call(rbind, res2)

Analisis discriminante con alta colinearidad

Datos originales

#### Disc datos originales

ad <- lda(sp ~ v1 + v2 + v3, data = res2T, CV = T)

tab1 <- table(res2T$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G  20   0   0   0   0   0
##   F1p   0  20   0   0   0   0
##   F2G   0   0  20   0   0   0
##   F2p   0   0   0  20   0   0
##   F3G   0   0   0   0  20   0
##   F3p   0   0   0   0   0  20
diag(prop.table(tab1,1))
## F1G F1p F2G F2p F3G F3p 
##   1   1   1   1   1   1
sum(diag(prop.table(tab1)))
## [1] 1
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(res2T$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
## F1 F2 F3 
##  1  1  1
sum(diag(prop.table(tab2)))
## [1] 1
cor(res2T[,-1])
##           v1        v2        v3
## v1 1.0000000 0.7059680 0.5482237
## v2 0.7059680 1.0000000 0.2022557
## v3 0.5482237 0.2022557 1.0000000

proporcionales al promedio de cada variable dentro de especies

#### Disc proporcionales al promedio de cada variable dentro de especies
dpvs <- lapply(unique(res2T$sp), function(x){
y <- res2T[res2T$sp == x, ] 

y[,-1] <- y[,-1]/colMeans(y[,-1])
return(y)
})

dpvs <- do.call(rbind, dpvs)


ad <- lda(sp ~ v1 + v2 + v3, data = dpvs, CV = T)

tab1 <- table(dpvs$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G   0  18   0   2   0   0
##   F1p  16   0   0   2   2   0
##   F2G   7   2   0   6   5   0
##   F2p   7   3   5   0   3   2
##   F3G   0   6   0   0   0  14
##   F3p   0   6   0   0  14   0
diag(prop.table(tab1,1))
## F1G F1p F2G F2p F3G F3p 
##   0   0   0   0   0   0
sum(diag(prop.table(tab1)))
## [1] 0
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dpvs$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
##    F1    F2    F3 
## 0.850 0.275 0.700
sum(diag(prop.table(tab2)))
## [1] 0.6083333

proporcionales a una medida dentro de especie

#### Disc proporcionales a una medida dentro de especie
#prop a v1
dv1s <- lapply(unique(res2T$sp), function(x){
  y <- res2T[res2T$sp == x, ] 
  
  y[,-1] <- y[,-1]/y$v2
  return(y)
})

dv1s <- do.call(rbind, dv1s)


ad <- lda(sp ~  v1 + v3, data = dv1s, CV = T)

tab1 <- table(dv1s$sp, ad$class)

tab1
##      
##       F1G F1p F2G F2p F3G F3p
##   F1G  13   7   0   0   0   0
##   F1p  12   8   0   0   0   0
##   F2G   0   0  13   7   0   0
##   F2p   0   0  16   4   0   0
##   F3G   0   0   0   0  15   5
##   F3p   0   0   0   0  10  10
diag(prop.table(tab1,1))
##  F1G  F1p  F2G  F2p  F3G  F3p 
## 0.65 0.40 0.65 0.20 0.75 0.50
sum(diag(prop.table(tab1)))
## [1] 0.525
cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dv1s$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
## F1 F2 F3 
##  1  1  1
sum(diag(prop.table(tab2)))
## [1] 1
cor(dv1s[,-c(1, 3)])
##           v1        v3
## v1 1.0000000 0.4853474
## v3 0.4853474 1.0000000

proporcionales a tamano de cada especie

#### 

dTSs <- lapply(unique(res2T$sp), function(x){
  y <- res2T[res2T$sp == x, ] 
  y[,-1] <- y[,-1]/mean(unlist(y[,-1]))
  return(y)
})

dTSs <- do.call(rbind, dTSs)


ad <- lda(sp ~  v1 + v2 + v3, data = dTSs, CV = T)

# tab1 <- table(dTSs$sp, ad$class)
# 
# tab1
# 
# diag(prop.table(tab1,1))
# 
# sum(diag(prop.table(tab1)))
# 

cat("clasificacion al nivel de forma")
## clasificacion al nivel de forma
tab2 <- table(substring(dTSs$sp, 1, 2), substring(ad$class, 1, 2))

diag(prop.table(tab2,1))
## F1 F2 F3 
##  1  1  1
sum(diag(prop.table(tab2)))
## [1] 1
cor(dTSs[,-1])
##           v1         v2         v3
## v1 1.0000000  0.1457167  0.2198478
## v2 0.1457167  1.0000000 -0.6044368
## v3 0.2198478 -0.6044368  1.0000000

Datos reales

proporcionales a una medida dentro de especie

nidos <- read.table(file = "/home/ba/Downloads/nidos.txt", header = T)

nidoprop <- lapply(unique(nidos$sp), function(x){
  y <- nidos[nidos$sp == x, ] 
  
  y[,-1] <- y[,-1]/y$vel
  return(y[,-grep("vel$",colnames(y))])
})

nidoprop <- do.call(rbind, nidoprop)

library(corrplot)


corrplot(cor(nidos[,-1]), main = "datos originales", mar = rep(3,4))

mean(cor(nidos[,-1])[upper.tri(cor(nidos[,-1]))])
## [1] 0.6496163
corrplot(cor(nidoprop[,-1]), main = " proporcionales a vel", mar = rep(3,4))

mean(cor(nidoprop[,-1])[upper.tri(cor(nidoprop[,-1]))])
## [1] 0.3782449

proporcionales a tamano promedio de nido

nidoprop2 <-  lapply(unique(nidos$sp), function(x){
  y <- nidos[nidos$sp == x, ] 
  y[,-1] <- y[,-1]/mean(unlist(y[,-1]))
  return(y)
})

nidoprop2 <- do.call(rbind, nidoprop2)

corrplot(cor(nidoprop2[,-1]), main = "proporcionales a tamano promedio de nido", mar = rep(3,4))

mean(cor(nidoprop2[,-1])[upper.tri(cor(nidoprop2[,-1]))])
## [1] 0.1902228

Por especie

for(i in unique(nidoprop2$sp))
  {
corrplot(cor(nidoprop2[nidoprop2$sp == i,-1]), main = i, mar = rep(3,4))
mean(cor(nidoprop2[nidoprop2$sp == i,-1])[upper.tri(cor(nidoprop2[nidoprop2$sp == i,-1]))])
  
  }