Fluorometro y Clorometro

Datos

primero se cargan los datos

library(readxl)
lab2 <- read_excel("C:/Users/Javi3/Downloads/lab2.xlsx")

Luego se crean los Vectores a trabajar

exf <- lab2$F[1:24] # Datos de Fv/m en el exterior 
inf <- lab2$F[25:48] # Datos del Fv/m en el interior
exc <- lab2$Clorofila[1:24] # Datos de clorofila en el exterior 
inc <- lab2$Clorofila[25:48] # Datos del clorofila en el interior

Supuestos

Normalidad

Ahora se determina por medio de la prueba Shapiro si hay Normalidad en los datos

shapiro.test(exc) #Test shapiro-wilks para los datos de Clorofila en el exterior

## 
##  Shapiro-Wilk normality test
## 
## data:  exc
## W = 0.94368, p-value = 0.1969

shapiro.test(inc) #Test shapiro-wilks para los datos de Clorofila en el interior

## 
##  Shapiro-Wilk normality test
## 
## data:  inc
## W = 0.97258, p-value = 0.7306

shapiro.test(exf) #Test shapiro-wilks para los datos de Fv/m en el exterior

## 
##  Shapiro-Wilk normality test
## 
## data:  exf
## W = 0.96415, p-value = 0.5273

shapiro.test(inf) #Test shapiro-wilks para los datos de Fv/m en el interior

## 
##  Shapiro-Wilk normality test
## 
## data:  inf
## W = 0.43154, p-value = 1.3e-08

Residuos

Pero lo que se debe evaluar es si los residuos siquen una distribucion normal por eso se hace una ligera modificacion

an <- lm(Clorofila ~ Ubi, data = lab2) #Anova para calcular los residuos de clorofila 
shapiro.test(residuals(an)) #Prueba Shapiro a los residuos de la clorofila

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(an)
## W = 0.9628, p-value = 0.131

an1 <- lm(F ~ Ubi, data = lab2) #Anova para calcular los residuos de Fv/m
shapiro.test(residuals(an1)) #Prueba Shapiro para los residuos de F

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(an1)
## W = 0.88846, p-value = 0.0002722

se hacen qqplots para analizar la normalidad de forma Grafica

qqnorm(resid(an)) # qqplot de la clorofila
qqline(resid(an))

qqnorm(resid(an1)) #qqplot de fv/m
qqline(resid(an1))

Homocedasticidad

ahora se analiza la Homogeneidad de las muestras evaluadad

var.test(exc,inc) #Prueba de homocedasticida para la clorofila

## 
##  F test to compare two variances
## 
## data:  exc and inc
## F = 1.2142, num df = 23, denom df = 23, p-value = 0.6454
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5252752 2.8069036
## sample estimates:
## ratio of variances 
##           1.214247

como los datos de Fv/m no cumplen con la normalidad se hace la prueba de Levene y la de Brown-Forsy

library(car)
leveneTest(lab2$F ~ as.factor(lab2$Ubi),center = mean) #Prueba de homocedasticidad para Fv/m

## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value    Pr(>F)    
## group  1  24.132 1.179e-05 ***
##       46                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

leveneTest(lab2$F ~ as.factor(lab2$Ubi), center = median) #Prueba de Brown-Forsy para Fv/m

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value    Pr(>F)    
## group  1  24.248 1.133e-05 ***
##       46                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

para analizar la homocedasticidad de forma Grafica

   ec <- resid(an)
   ef <- resid(an1)
   rec<-rstandard(an)
   ref<-rstandard(an1)
   prec<-predict(an)
   pref<-predict(an1)
   plot(prec, rec, xlab="Predichos de Clorofila", ylab="Residuos estandarizados de Clorofila",main="Gráfico de dispersión de RE vs PRED de Clorofila" )
   abline(0,0)

   plot(pref, ref, xlab="Predichos Fv/m", ylab="Residuos estandarizados Fv/m",main="Gráfico de dispersión de RE vs PRED Fv/m" )
   abline(0,0)

Analisis

Para los datos de la clorofila que si cumple con los supuestos, por lo que es posible realizar una prueba t-Student

t.test(exc, inc, var.equal = TRUE) #Prueba para analizar las diferencias entre la clorofila

## 
##  Two Sample t-test
## 
## data:  exc and inc
## t = -5.9437, df = 46, p-value = 3.516e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -11.776593  -5.818036
## sample estimates:
## mean of x mean of y 
##  48.94296  57.74028

para analizar los datos del Fvm se realiza un Bootstrap no parametrico

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:car':
## 
##     logit

set.seed(0728)
x <- exf - inf   #Simulacion de la diferencia
sampleMean <- rep(NA, 15000) # vectores para que el boot 
sampleVar <- rep(NA, 15000)
for(i in 1:15000) {sampleMean [i] <-mean(rnorm(length(x), mean(x), sd(x))) # pone valores a los vectores  
sampleVar [i] <- var(rnorm(length(x),mean(x), sd(x)))}
mediaFunc <- function(x, i) { return(c(mean(x[i]), var(x[i])))} # La funcion del boot  
bootMedia <- boot(x, mediaFunc, 15000) # boot con 15000 replicas
summary(bootMedia)

##       R  original    bootBias    bootSE   bootMed
## 1 15000 -0.378671  0.00027904 0.0375071 -0.378921
## 2 15000  0.034826 -0.00144320 0.0087771  0.032952

plot(bootMedia) # plot del boot

boot.ci(bootMedia) # Intervalos de confianza del boot

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 15000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bootMedia)
## 
## Intervals : 
## Level      Normal              Basic             Studentized     
## 95%   (-0.4525, -0.3054 )   (-0.4534, -0.3057 )   (-0.4547, -0.2952 )  
## 
## Level     Percentile            BCa          
## 95%   (-0.4516, -0.3039 )   (-0.4490, -0.3018 )  
## Calculations and Intervals on Original Scale

igual mente se realiza una simulacion Bayesiana con el mcmc

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.2.0

## Loaded modules: basemod,bugs

set.seed(0728)
sef <- rnorm(length(exf), mean(exf), sd(exf))
sif <- rnorm(length(inf), mean(inf), sd(inf))
mcf <- sef - sif
#Especificaciones del modelo para el muestreo gibs
model_string <- "model{     
   for(i in 1:length(y)) {
       y[i] ~ dnorm(mu, tau)
   }
   mu ~ dnorm(0, tau/w)
   w <- tau/tau0
   invsigma ~ dgamma(0.000001, 0.000001)
   tau <- 1 / invsigma
   tau0 <- 90
 }"



#Para correr el modelo con jags
modelo <- jags.model(textConnection(model_string), data = list(y = mcf), n.chains = 3, n.adapt = 15000)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 24
##    Unobserved stochastic nodes: 2
##    Total graph size: 36
## 
## Initializing model

update(modelo,15000)
mcmc_muestra <- coda.samples(modelo,variable.names= c("mu", "invsigma"), n.iter = 30000)
plot(mcmc_muestra) # Muestra las densidades marginales

summary(mcmc_muestra) #Para obtener los percentiles bayesianos

## 
## Iterations = 30001:60000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 30000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##              Mean      SD  Naive SE Time-series SE
## invsigma  0.04331 0.01555 5.182e-05      0.0001040
## mu       -0.36619 0.04270 1.423e-04      0.0001865
## 
## 2. Quantiles for each variable:
## 
##              2.5%      25%      50%      75%    97.5%
## invsigma  0.02289  0.03266  0.04012  0.05033  0.08205
## mu       -0.44339 -0.39542 -0.36872 -0.33954 -0.27576