#CARPETA ORIGINAL:
#  C:\Users\Majo\Documents\POSDOC\A_Manuscripto\R\C. pomonella: 
#Analisis de proporcion
#Del R 
#FECHA DE TRABAJO OFICIAL 
#Analisis de proporciones de mortalidad entre densidad.
#DATOS actuales en: c:\Users\Majo\Documents\POSDOC\A_Manuscripto\R\C. pomonella\Analisis de proporcion
#TABLA DE DATOS: "Tabla General DATOS_ CP-GL_CMJ5"
#TABA CSV: "datos"
# Comparacion de  proporciones
# Se analiza sobre el total de densidad evaluadas.
# Test de igualdad de proporciones. Funci?n: prop.test ().
# Conclusion: con un error de tipo I menor a 0,05 concluimos 
# que las proporciones de mortalirdad entre  grupos varia.
# cuales son los grupos que difieren. Funcion pairwise.prop.test ()

rm(list=ls())
ls()

#Carga de datos


tablaMortalidad<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
tablaPARASITISMO<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
tablaHF<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")


tablaMortalidad
tablaPARASITISMO
tablaHF


> tablaMortalidad

1     8      8
2    16     15
3    32     19
4    64     32
5   128     52


#Calculo de proporciones por  densidad

prop.test(tablaMortalidad$efecto,tablaMortalidad$total,conf.level=0.95)


prop.test(tablaMortalidad$efecto,tablaMortalidad$total,conf.level=0.95)

5-sample test for equality of proportions without continuity correction

data:  tablaMortalidad$efecto out of tablaMortalidad$total
X-squared = 25.817, df = 4, p-value = 3.445e-05
alternative hypothesis: two.sided
sample estimates:
  prop 1  prop 2  prop 3  prop 4  prop 5 
1.00000 0.93750 0.59375 0.50000 0.40625  

#Analisis estadistico

pairwise.prop.test(tablaMortalidad$efecto,tablaMortalidad$total,conf.level=0.95)



Pairwise comparisons using Pairwise comparison of proportions 

data:  tablaMortalidad$efecto out of tablaMortalidad$total 

1      2      3      4     
2 1.0000 -      -      -     
  3 0.3817 0.1975 -      -     
  4 0.1478 0.0321 1.0000 -     
  5 0.0321 0.0018 0.3817 0.8403



#PARASITISMO
tablaPARASITISMO
prop.test(tablaPARASITISMO$Parasitismo,tablaPARASITISMO$Total,conf.level=0.95)
pairwise.prop.test(tablaPARASITISMO$Parasitismo,tablaPARASITISMO$Total,conf.level=0.95)

#HOST FEEDING
tablaHF
tablaHF<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
tablaHF
prop.test(tablaHF$HF,tablaHF$Total,conf.level=0.95)
pairwise.prop.test(tablaHF$HF,tablaHF$Total,conf.level=0.95)


#Comparaci?n por densidad, entre parasitismo y HF
#ANALISIS DE OR


# NO MODIFICAR
# aqui comienza la definicion de la funcion 
# calcOddsRatio. #el calculo se hace con el dato
# de la fila 1 como referencia (referencerow)
calcOddsRatio <- function(mymatrix,alpha=0.05,referencerow=2,quiet=FALSE)
{
  numrow <- nrow(mymatrix)
  myrownames <- rownames(mymatrix)
  
  for (i in 1:numrow)
  {
    rowname <- myrownames[i]
    DiseaseUnexposed <- mymatrix[referencerow,1]
    ControlUnexposed <- mymatrix[referencerow,2]
    if (i != referencerow)
    {
      DiseaseExposed <- mymatrix[i,1]
      ControlExposed <- mymatrix[i,2]
      
      totExposed <- DiseaseExposed + ControlExposed
      totUnexposed <- DiseaseUnexposed + ControlUnexposed
      
      probDiseaseGivenExposed <- DiseaseExposed/totExposed
      probDiseaseGivenUnexposed <- DiseaseUnexposed/totUnexposed
      probControlGivenExposed <- ControlExposed/totExposed
      probControlGivenUnexposed <- ControlUnexposed/totUnexposed
      
      # calculate the odds ratio
      oddsRatio <- (probDiseaseGivenExposed*probControlGivenUnexposed)/
        (probControlGivenExposed*probDiseaseGivenUnexposed)
      if (quiet == FALSE)
      {
        print(paste("category =", rowname, ", odds ratio = ",oddsRatio))
      }
      
      # calculate a confidence interval
      confidenceLevel <- (1 - alpha)*100
      sigma <- sqrt((1/DiseaseExposed)+(1/ControlExposed)+
                      (1/DiseaseUnexposed)+(1/ControlUnexposed))
      # sigma is the standard error of our estimate of the log of the odds ratio
      z <- qnorm(1-(alpha/2))
      lowervalue <- oddsRatio * exp(-z * sigma)
      uppervalue <- oddsRatio * exp( z * sigma)
      if (quiet == FALSE)
      {
        print(paste("category =", rowname, ", ", confidenceLevel,
                    "% confidence interval = [",lowervalue,",",uppervalue,"]"))
      }
    }
  }
  if (quiet == TRUE && numrow == 2) # If there are just two treatments (exposed/nonexposed)
  {
    return(oddsRatio)
  }
}
# FIN definicion de funcion


# Datos: c:\Users\Majo\Documents\POSDOC\2020\A_Manuscripto\R\C. pomonella\Analisis de proporcion
#solapa:OR tabla R

# comparaci?n odds ratio 
#mymatrix tiene la matriz de datos a analizar
#mymatrix <- matrix(c(A,B,C,D),nrow=2,byrow=TRUE)
#A=  + en PARASITISMO
#B=  - en  PARASITISMO
#C=  + en HF
#D= - en HF
#Si el intervalo incluye el 1: NO SIGNIFICATIVO

# DENSIDAD 1

# Comparacion de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(6,2,2,6),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 2

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(12,4,3,13),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

#Densidad 4
# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(12,20,7,25),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 8

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(13,51,19,45),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 16

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(13,115,39,89),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

####
#Comparacion de HF con parasitismo / MISMOS RESULTADOS
#Comparaci?n por densidad, entre parasitismo y HF
#ANALISIS DE OR
# NO MODIFICAR
# aqui comienza la definicion de la funcion 
# calcOddsRatio. #el calculo se hace con el dato
# de la fila 1 como referencia (referencerow)
calcOddsRatio <- function(mymatrix,alpha=0.05,referencerow=2,quiet=FALSE)
{
  numrow <- nrow(mymatrix)
  myrownames <- rownames(mymatrix)
  
  for (i in 1:numrow)
  {
    rowname <- myrownames[i]
    DiseaseUnexposed <- mymatrix[referencerow,1]
    ControlUnexposed <- mymatrix[referencerow,2]
    if (i != referencerow)
    {
      DiseaseExposed <- mymatrix[i,1]
      ControlExposed <- mymatrix[i,2]
      
      totExposed <- DiseaseExposed + ControlExposed
      totUnexposed <- DiseaseUnexposed + ControlUnexposed
      
      probDiseaseGivenExposed <- DiseaseExposed/totExposed
      probDiseaseGivenUnexposed <- DiseaseUnexposed/totUnexposed
      probControlGivenExposed <- ControlExposed/totExposed
      probControlGivenUnexposed <- ControlUnexposed/totUnexposed
      
      # calculate the odds ratio
      oddsRatio <- (probDiseaseGivenExposed*probControlGivenUnexposed)/
        (probControlGivenExposed*probDiseaseGivenUnexposed)
      if (quiet == FALSE)
      {
        print(paste("category =", rowname, ", odds ratio = ",oddsRatio))
      }
      
      # calculate a confidence interval
      confidenceLevel <- (1 - alpha)*100
      sigma <- sqrt((1/DiseaseExposed)+(1/ControlExposed)+
                      (1/DiseaseUnexposed)+(1/ControlUnexposed))
      # sigma is the standard error of our estimate of the log of the odds ratio
      z <- qnorm(1-(alpha/2))
      lowervalue <- oddsRatio * exp(-z * sigma)
      uppervalue <- oddsRatio * exp( z * sigma)
      if (quiet == FALSE)
      {
        print(paste("category =", rowname, ", ", confidenceLevel,
                    "% confidence interval = [",lowervalue,",",uppervalue,"]"))
      }
    }
  }
  if (quiet == TRUE && numrow == 2) # If there are just two treatments (exposed/nonexposed)
  {
    return(oddsRatio)
  }
}
# FIN definicion de funcion


# Datos: c:\Users\Majo\Documents\POSDOC\2020\A_Manuscripto\R\C. pomonella\Analisis de proporcion
#solapa:OR tabla R

# comparacion odds ratio 
#mymatrix tiene la matriz de datos a analizar
#mymatrix <- matrix(c(A,B,C,D),nrow=2,byrow=TRUE)
#A=  + en PARASITISMO
#B=  - en  PARASITISMO
#C=  + en HF
#D= - en HF
#Si el intervalo incluye el 1: NO SIGNIFICATIVO

# DENSIDAD 1

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(2,6,6,2),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 2

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(3,13,12,4),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 4

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(7,25,12,20),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 8

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(19,45,13,51),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)

# DENSIDAD 16

# Comparaci?n de datos
# negparat,positivparasit, negHF, positHF

mymatrix <- matrix(c(39,89,13,115),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)



mymatrix <- matrix(c(6,2,2,6),nrow=2,byrow=TRUE) 
colnames(mymatrix) <- c("Pos","Neg")
rownames(mymatrix) <- c("Parasitismo","HF")
print(mymatrix)
calcOddsRatio(mymatrix,alpha=0.05)



##Datos-Respuesta funcional-mean y desvio
#Tabla 1
#25/04/2023
#Sumar a tabla 1  medias y desvio estandar  de valores de 
#mortalidad por densidad
#Documents\Carrera\2023\ms 2023_ R. funcional\
#Estadictos_Nuevo analisis
# Tabla de datos 
# Vector de muestra
#Cada vector es la suma de mortalidad por densidad por replica

Dens1 <- c(1, 1, 1, 1, 1, 1, 1, 1)
Dens2 <- c(1, 2, 2, 2, 2, 2, 2, 2)
Dens4 <- c(4, 2, 2, 2, 2, 4, 1, 2)
Dens8 <- c(6, 4, 5, 3, 6, 6, 1, 1)
Dens16 <- c(7, 6, 11, 9, 6, 4, 2,7)

# Desviacion tipica
sd(Dens1)
sd(Dens2)
sd(Dens4)
sd(Dens8)
sd(Dens16)


mean(Dens1)
mean(Dens2)
mean(Dens4)
mean(Dens8)
mean(Dens16)    



#Parasitismo
# Vector de muestra: Parasitismo
pDens1 <- c(1, 0, 1, 1, 1, 0, 1, 1)
pDens2 <- c(2, 1, 1, 1, 2, 2, 2, 1)
pDens4 <- c(4, 1, 0, 0, 2, 2, 1, 2)
pDens8 <- c(4, 2, 3, 0, 1, 2, 1, 0)
pDens16 <- c(2, 1, 4, 2, 1, 2, 1, 0)

# Desviaci?n t?pica
sd(pDens1)
sd(pDens2)
sd(pDens4)
sd(pDens8)
sd(pDens16)

#Media
mean(pDens1)
mean(pDens2)
mean(pDens4)
mean(pDens8)
mean(pDens16)


#Response functional
#Tabla Base: RF 

install.packages("frair")               
library("frair")

#Table: RF
RF<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
View(RF)
names(RF)
frair_test(formula = efecto~total, data = Base)
frair_test(formula = muertos~densidad, data = RF)


#Interferencia mutua
#CARPETA ORIGINAL:
#  C:\Users\Majo\Documents\POSDOC\2020\A_Manuscripto\R\C. pomonella: 
#Analisis de cantidad de huevos: individualizacion de larvas
#Del R: 
#FECHA DE TRABAJO OFICIAL 
#Analisis 
#DATOS actuales en: c:\Users\Majo\Documents\carrerra\A_Manuscripto\R\C. pomonella\Analisis de proporcion
#TABLA DE DATOS: "Tabla General DATOS_ CP-GL_CMJ5"
#TABA CSV: "datos"
# Se analiza sobre el total de densidad evaluadas.
#anava
#Datos de progenie
#Nombre de tabla: ProgenieIM

ProgenieIM<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")

names(ProgenieIM)

#Anova por Tratamiento
#Anova "Individualizaci?n.de.huevos.por.larva."

anovaHuevos <- aov(Huevos ~ Tratamiento, ProgenieIM)
summary(anovaHuevos)
TukeyHSD(anovaHuevos)

anovaAdultos <- aov(Adultos ~ Tratamiento, ProgenieIM)
summary(anovaAdultos)
TukeyHSD(anovaAdultos)

anovaRelacionHyA<- aov(RelacionHyA ~ Tratamiento, ProgenieIM)
summary(anovaRelacionHyA)
TukeyHSD(anovaRelacionHyA)


anovaProporcionEmergencia <- aov(ProporcionEmergencia  ~ Tratamiento, ProgenieIM)
summary(anovaProporcionEmergencia)
TukeyHSD(anovaProporcionEmergencia)


install.packages("ggplot2")
library(ggplot2)


ggplot(data = ProgenieIM, aes(x = Huevos, y = Tratamiento, color = Tratamiento)) +
  geom_boxplot() +
  theme_bw()


ggplot(data = ProgenieIM, aes(x = Adultos, y = Tratamiento, color = Tratamiento)) +
  geom_boxplot() +
  theme_bw()

ggplot(data = ProgenieIM, aes(x = RelacionHyA, y = Tratamiento, color = Tratamiento)) +
  geom_boxplot() +
  theme_bw()


ggplot(data = ProgenieIM, aes(x = ProporcionEmergencia , y = Tratamiento, color = Tratamiento)) +
  geom_boxplot() +
  theme_bw()


#DESVIO ESTANDAR
#Tabla: datos de progenie
#Tabla: sd, datos de progenie
#Tabla: ProgenieIM
#Datos larvas parasitadas
#Por tratamiento

#Datos 
#Tratamiento A
#Tabla: ProgenieIM-A
#Tabla: ProgenieIM-B
#Tabla: ProgenieIM-C



ProgenieIMA<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
ProgenieIMB<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")
ProgenieIMC<-read.table("clipboard",header=TRUE,dec=",",sep="\t",encoding="latin1")


#Tratamiento A
sd(ProgenieIMA$Huevos)
sd(ProgenieIMA$Adultos)
sd(ProgenieIMA$RelacionHyA)
sd(ProgenieIMA$ProporcionEmergencia)

#Tratamiento B
sd(ProgenieIMB$Huevos)
sd(ProgenieIMB$Adultos)
sd(ProgenieIMB$RelacionHyA)
sd(ProgenieIMB$ProporcionEmergencia )


#Tratamiento C
sd(ProgenieIMC$Huevos)
sd(ProgenieIMC$Adultos)
sd(ProgenieIMC$RelacionHyA)
sd(ProgenieIMC$ProporcionEmergencia)


#fIN del script
#No modificar