IRGA

Datos

library(readxl)
jav <- read_excel("C:/Users/Javi3/Downloads/irga.xlsx")
expa <- jav$PAR[1:24] # Datos de PAR en el exterior 
inpa <- jav$PAR[25:48] # Datos del PAR en el interior
exp <- jav$Pn[1:24] # Datos de Pn en el exterior
inp <- jav$Pn[25:48] # Datos del Pn en el interior
exe <- jav$E[1:24] # Datos de E en el exterior
ine <- jav$E[25:48] # Datos del E en el interior
exc <- jav$C[1:24] # Datos de C en el exterior
inc <- jav$C[25:48] # Datos del C en el interior

Normalidad

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

## 
##  Shapiro-Wilk normality test
## 
## data:  expa
## W = 0.92022, p-value = 0.05906

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

## 
##  Shapiro-Wilk normality test
## 
## data:  inpa
## W = 0.34702, p-value = 2.471e-09

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

## 
##  Shapiro-Wilk normality test
## 
## data:  exp
## W = 0.5887, p-value = 4.631e-07

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

## 
##  Shapiro-Wilk normality test
## 
## data:  inp
## W = 0.64649, p-value = 2.134e-06

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

## 
##  Shapiro-Wilk normality test
## 
## data:  exe
## W = 0.95272, p-value = 0.31

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

## 
##  Shapiro-Wilk normality test
## 
## data:  ine
## W = 0.97286, p-value = 0.7376

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

## 
##  Shapiro-Wilk normality test
## 
## data:  exc
## W = 0.95676, p-value = 0.3768

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

## 
##  Shapiro-Wilk normality test
## 
## data:  inc
## W = 0.47488, p-value = 3.244e-08

Homocedasticidad

var.test(exe,ine) # Prueba de Homocedasticidad para muestras que cumplen el supuesto de normalidad

## 
##  F test to compare two variances
## 
## data:  exe and ine
## F = 3.217, num df = 23, denom df = 23, p-value = 0.006953
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.391674 7.436661
## sample estimates:
## ratio of variances 
##           3.217049

library(car)
leveneTest(jav$PAR ~ as.factor(jav$Ubi)) #Prueba de homocedasticidad para PAR

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

leveneTest(jav$Pn ~ as.factor(jav$Ubi)) #Prueba de homocedasticidad para Pn

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1494 0.7009
##       46

leveneTest(jav$C ~ as.factor(jav$Ubi)) #Prueba de homocedasticidad para C

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

Pruebas Estadistica

Parametrica

t.test(exe,ine) #Prueba t student para dos muestras con varianzas diferentes

## 
##  Welch Two Sample t-test
## 
## data:  exe and ine
## t = 0.011588, df = 36.039, p-value = 0.9908
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1933503  0.1955726
## sample estimates:
## mean of x mean of y 
##  1.204306  1.203194

No Parametrica

wilcox.test(exp, inp) #Prueba Wilcox para Pn

## 
##  Wilcoxon rank sum test
## 
## data:  exp and inp
## W = 434, p-value = 0.002184
## alternative hypothesis: true location shift is not equal to 0

No Deterministas

Bootstrap

remuestreo para C

library(boot)

## 
## Attaching package: 'boot'

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

x1 <- exc - inc 
sampleMean <- rep(NA, 10000)
for(i in 1:10000) {sampleMean [i] <- mean(x1)}
meanFunc <- function(x1, i) { mean(x1[i])}
bootMean <- boot(x1, meanFunc, 10000)
plot(bootMean)

summary(bootMean)

##       R original bootBias bootSE bootMed
## 1 10000  -1333.8   0.4146 554.65 -1291.2

boot.ci(bootMean)

## Warning in boot.ci(bootMean): bootstrap variances needed for studentized
## intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bootMean)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-2421,  -247 )   (-2302,  -131 )  
## 
## Level     Percentile            BCa          
## 95%   (-2537,  -366 )   (-2938,  -555 )  
## Calculations and Intervals on Original Scale

Remuestreo para PAR

x <- expa - inpa 
sampleMean1 <- rep(NA, 10000)
for(i in 1:10000) {sampleMean [i] <- mean(x)}
bootMean1 <- boot(x, meanFunc, 10000)
plot(bootMean1)

summary(bootMean1)

##       R original bootBias bootSE bootMed
## 1 10000   652.38 -0.13653 104.41  650.68

boot.ci(bootMean1)

## Warning in boot.ci(bootMean1): bootstrap variances needed for studentized
## intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bootMean1)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (447.9, 857.2 )   (447.0, 854.0 )  
## 
## Level     Percentile            BCa          
## 95%   (450.8, 857.8 )   (456.4, 863.7 )  
## Calculations and Intervals on Original Scale

Monte Carlo

Simulacion de la diferencia entre las medias de C

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.2.0

## Loaded modules: basemod,bugs

set.seed(95)
sec <- rt(100, 47, mean(exc))
sic <- rt(100, 47, mean(inc))
mcc <- sec - sic
model_string <- "model{
   for(i in 1:length(y)) {
     y[i] ~ dnorm(mu, tau)
   }
   mu ~ dnorm(0, 0.0001)
   sigma ~ dlnorm(0, 0.0625)
   tau <- 1 / pow(sigma, 2)
 }"

model <- jags.model(textConnection(model_string), data = list(y = mcc), n.chains = 3, n.adapt = 10000)

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

update(model, 10000)
mcmc_samples <- coda.samples(model, variable.names=c("mu", "sigma"), n.iter = 20000)
plot(mcmc_samples)

summary(mcmc_samples)

## 
## Iterations = 20001:40000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 20000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean    SD Naive SE Time-series SE
## mu    -1350.1 15.37  0.06274        0.07067
## sigma   148.4 11.11  0.04534        0.06395
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%     75%   97.5%
## mu    -1379.2 -1360.6 -1350.4 -1340.1 -1319.0
## sigma   128.8   140.6   147.7   155.4   172.2

Simulacion de la diferencia para PAR

sepa <- rt(100, 47, mean(expa))
sipa <- rt(100, 47, mean(inpa))
mcpa <- sepa - sipa
model_string1 <- "model{
   for(i in 1:length(y)) {
     y[i] ~ dnorm(mu, tau)
   }
   mu ~ dnorm(0, 0.0001)
   sigma ~ dlnorm(0, 0.0625)
   tau <- 1 / pow(sigma, 2)
 }"

model <- jags.model(textConnection(model_string1), data = list(y = mcpa), n.chains = 3, n.adapt = 10000)

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

update(model, 10000)
mcmc_samples1 <- coda.samples(model, variable.names=c("mu", "sigma"), n.iter = 20000)
plot(mcmc_samples1)

summary(mcmc_samples1)

## 
## Iterations = 20001:40000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 20000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean    SD Naive SE Time-series SE
## mu    648.26 6.727  0.02746        0.02775
## sigma  67.08 4.839  0.01976        0.02602
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%    50%    75%  97.5%
## mu    634.93 643.74 648.27 652.82 661.32
## sigma  58.41  63.73  66.81  70.14  77.38