Script MixSIAR - Ana Maria Kruger

---

DATA

library(MixSIAR)
getOption("max.print")

[1] 9999999

options (max.print = 9999999)#esta função simplesmente mostra o numero de elementos as ser mostrado na tela

# Load mixture data
mix.filename = system.file("extdata", "dataANA_consumer.csv", package = "MixSIAR")
mix = load_mix_data(filename=mix.filename,iso_names=c("d13c","d15N"),factors=c("id"),fac_random=c(FALSE), 
      fac_nested=c(FALSE),cont_effects=NULL)
mix

$data

$data_iso
        d13c  d15N
 [1,] -31.91  7.36
 [2,] -31.86  7.44
 [3,] -28.33  5.71
 [4,] -27.61  5.66
 [5,] -24.22  4.40
 [6,] -24.23  4.52
 [7,] -19.98  2.52
 [8,] -20.65  2.69
 [9,] -16.28  1.34
[10,] -16.37  0.93
[11,] -12.92 -0.73
[12,] -12.79 -0.52

$n.iso
[1] 2

$n.re
[1] 0

$n.ce
[1] 0

$FAC
$FAC[[1]]
$FAC[[1]]$values
 [1] 2 2 6 6 5 5 4 4 3 3 1 1

$FAC[[1]]$levels
[1] 6

$FAC[[1]]$labels
[1] "100B"   "100T"   "20T80B" "40T60B" "60T40B" "80T20B"

$FAC[[1]]$lookup
NULL

$FAC[[1]]$re
[1] FALSE

$FAC[[1]]$name
[1] "id"



$CE
list()

$CE_orig
list()

$CE_center
logical(0)

$CE_scale
logical(0)

$cont_effects
NULL

$MU_names
[1] "Meand13c" "Meand15N"

$SIG_names
[1] "SDd13c" "SDd15N"

$iso_names
[1] "d13c" "d15N"

$N
[1] 12

$n.fe
[1] 1

$n.effects
[1] 1

$factors
[1] "id"

$fac_random
[1] FALSE

$fac_nested
[1] FALSE

$fere
[1] FALSE

# Load source data
source.filename = system.file("extdata", "dataANA_sources.csv", package = "MixSIAR")
source = load_source_data(filename=source.filename, source_factors=NULL, conc_dep=FALSE, data_type="means", mix)
source

$n.sources
[1] 3

$source_names
[1] "Brachiaria" "Tithonia"   "Urea"      

$S_MU
  Meand13c Meand15N
1  -12.820    3.556
2  -30.994    8.553
3  -42.720   -0.625

$S_SIG
  SDd13c SDd15N
1  0.417  1.156
2  0.241  2.693
3  0.050  0.145

$S_factor1
NULL

$S_factor_levels
NULL

$conc
NULL

$MU_array
  Meand13c Meand15N
1  -12.820    3.556
2  -30.994    8.553
3  -42.720   -0.625

$SIG2_array
    SDd13c   SDd15N
1 0.173889 1.336336
2 0.058081 7.252249
3 0.002500 0.021025

$n_array
[1] 5 5 2

$SOURCE_array
NULL

$n.rep
NULL

$by_factor
[1] NA

$data_type
[1] "means"

$conc_dep
[1] FALSE

# Load discrimation data
discr.filename = system.file("extdata", "dataANA_discrimination.csv", package = "MixSIAR")
discr = load_discr_data(filename=discr.filename, mix)
discr

$mu
           Meand13c Meand15N
Brachiaria        0        0
Tithonia          0        0
Urea              0        0

$sig2
           SDd13c SDd15N
Brachiaria      0      0
Tithonia        0      0
Urea            0      0

OBS.: Estes 3 dataframes precisam ser preparados em .csv e salvados na pasta “extdata” do MixSIAR no disco local C: Na pasta “extdata” terá varios outros dataframes que são exemplos que já vem com o pacote, isso é bem comum no R, cuidado para não confundi-los com o seu. O caminho para encontrar a pasta “extdata” será algo assim: C:-library\4.2 O erros geralmente são devido ao fato de não ter salvado os dataframes no local correto, ou devido ao nome das colunas de cada dataframe, que devem ser semelhantes aos “iso_names” do primeiro dataframe.

Plot da dispersão dos dados

g=plot_data(filename="isospace_plot",
          plot_save_pdf=TRUE,
          plot_save_png=FALSE,
          mix,source,discr,
          return_obj=T)
g

#modify plot using ggplot2 commands

library(ggplot2)

g1=g+scale_x_continuous(name=expression(paste(delta^{13},"C (‰)")),breaks=seq(-12,-45,-3))+
  scale_y_continuous(name=expression(paste(delta^{15},"N (‰)")), breaks=seq(-4,10,2))+
  theme(plot.title = element_blank(),
  panel.grid.major = element_blank(), 
  panel.grid.minor = element_blank(), 
  axis.line = element_line(colour = "black"),
  panel.background = element_rect(fill = "transparent"),
  legend.background = element_rect(fill = "transparent",size=0.5, linetype="solid",colour ="black"),
  axis.title=element_text(size=20), legend.position= c(0.01, 0.9),axis.text=element_text(size=14),
  legend.text=element_text(size=10))+
  annotate("text", size=3.8, x=-42.5, y=11,label="Treatments:")
g1

library(cowplot)
#Save plot
save_plot("Plot da dispersão.pdf",g1, ncol = 1, nrow = 1) 

site usando está ferramenta caso surja alguma duvida: https://rdrr.io/cran/MixSIAR/man/plot_data.html Recomendo fazer as modificações no grafico pelo ggplot2 antes da importação

Modelo

output_JAGS(jags.1, mix, source, output_options)

################################################################################
# Gelman-Rubin Diagnostic
################################################################################

Generally the Gelman diagnostic should be < 1.05

Out of 46 variables: 0 > 1.01
                      0 > 1.05
                      0 > 1.1

The worst variables are:
              Point est. Upper C.I.
ilr.fac1[6,1]   1.001542   1.005244
loglik[7]       1.001066   1.002189
ilr.fac1[5,1]   1.000849   1.001076
p.fac1[1,1]     1.000838   1.003077
p.global[1]     1.000838   1.003077
p.fac1[3,1]     1.000762   1.002868
p.fac1[3,2]     1.000713   1.002360
p.fac1[1,3]     1.000682   1.001983
p.global[3]     1.000682   1.001983
deviance        1.000555   1.002466

################################################################################
# Geweke Diagnostic
################################################################################

The Geweke diagnostic is a standard z-score, so we'd expect 5% to be outside +/-1.96
Number of variables outside +/-1.96 in each chain (out of 46):

       Chain 1 Chain 2 Chain 3
Geweke       1       0       5

################################################################################
# Summary Statistics
################################################################################

DIC = 28.7372

                     Mean    SD  2.5%    5%   25%   50%   75%   95% 97.5%
p.100B.Brachiaria   0.902 0.033 0.824 0.840 0.883 0.907 0.926 0.946 0.951
p.100T.Brachiaria   0.095 0.037 0.031 0.039 0.069 0.093 0.119 0.160 0.173
p.20T80B.Brachiaria 0.861 0.038 0.785 0.798 0.837 0.861 0.884 0.925 0.942
p.40T60B.Brachiaria 0.694 0.048 0.593 0.607 0.662 0.699 0.727 0.762 0.777
p.60T40B.Brachiaria 0.477 0.058 0.372 0.385 0.436 0.477 0.516 0.577 0.595
p.80T20B.Brachiaria 0.279 0.055 0.175 0.184 0.240 0.282 0.319 0.366 0.382
p.100B.Tithonia     0.060 0.026 0.022 0.026 0.041 0.055 0.074 0.110 0.123
p.100T.Tithonia     0.684 0.095 0.490 0.520 0.623 0.689 0.752 0.830 0.850
p.20T80B.Tithonia   0.070 0.050 0.006 0.009 0.029 0.059 0.103 0.168 0.185
p.40T60B.Tithonia   0.149 0.102 0.013 0.019 0.065 0.127 0.219 0.349 0.376
p.60T40B.Tithonia   0.363 0.138 0.083 0.126 0.267 0.362 0.467 0.587 0.605
p.80T20B.Tithonia   0.549 0.138 0.301 0.336 0.450 0.540 0.648 0.787 0.806
p.100B.Urea         0.038 0.017 0.014 0.016 0.026 0.035 0.047 0.071 0.081
p.100T.Urea         0.220 0.058 0.117 0.130 0.179 0.217 0.259 0.321 0.342
p.20T80B.Urea       0.069 0.036 0.007 0.011 0.040 0.070 0.096 0.127 0.137
p.40T60B.Urea       0.157 0.065 0.019 0.033 0.113 0.168 0.209 0.246 0.257
p.60T40B.Urea       0.159 0.084 0.014 0.024 0.095 0.160 0.219 0.303 0.327
p.80T20B.Urea       0.172 0.084 0.015 0.026 0.109 0.177 0.233 0.303 0.323
NULL

Além dos dataframes, é necessario salvar um arquivo em .txt na pasta “extdata”, o nome do arquivo é “MixSIAR_model.txt”, na verdade este é o modelo utilizado. Quem nunca rodou o MixSIAR terpa alguns erros no script neste setor, erros relacionados ao “rjags” package, isso é devido ao fato de estar faltando o software JAGS no seu computador, um software estatístico que é usado para realizar inferência estatística, usado para realizar a inferência bayesiana, que envolve a atualização das estimativas das quantidades desconhecidas com base em dados observados e conhecimento prévio. Este software pode ser baixado no seguinte link: http://www.sourceforge.net/projects/mcmc-jags/files