The present analysis focused on determining the density of mabuyas (Trachylepis atlantica) in different habitats on the island of Fernando de Noronha. The R package used for the present analysis is the RMark and the code is hosted in github (https://github.com/tati-micheletti/mabuyas). This .rmd file will guide you on how to reproduce the analysis. Make sure you also have program Mark and RStudio installed in your computer. Mark needs to be installed in the ProgramFiles folder (C:\ProgramFiles\).
clone or download, and click on the download ZIP to download the folder.mabuyas.Rproj in RStudio.Ctrl + Alt + R or clicking in Run and Run All in the tab above.RMark package in R
Analysis were performed using the Poisson-log normal mark-resight model.
## Warning: package 'RMark' was built under R version 3.4.4## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip## Warning: package 'xlsx' was built under R version 3.4.4## Loading required package: rJava## Warning: package 'rJava' was built under R version 3.4.4## Loading required package: xlsxjars## Warning: package 'xlsxjars' was built under R version 3.4.4If you don’t have the libraries above, install and run them again:
install.packages("RMark")
install.packages("xlsx")
install.packages("reproducible")First of all, to analyse your data, you need to bring it to R.
excel to be analysed by Mark:directory <- getwd()
ch <- read.xlsx(file = file.path(getwd(), "data/dataMabuya.xlsx"),
header = TRUE,
sheetName = "ch",
as.data.frame = TRUE)
ch$ch <- as.character(ch$ch) # Convert factor to character
covars <- read.xlsx(file = file.path(getwd(), "data/dataMabuya.xlsx"),
header = TRUE,
sheetName = "covars",
as.data.frame = TRUE)nocc <- 6 # Number of occasions
groupsNames <- unique(covars$colonies) # Names of groups
nGroups <- length(groupsNames) # Number of groups
unmarkedSeen <- matrix(covars[, "unmarkedSeen"],
nrow = nGroups,
ncol = nocc,
byrow = TRUE) # counts of unmarked individuals seen
markedUnidentified <- matrix(covars[, "markedUnidentified"],
nrow = nGroups,
ncol = nocc,
byrow = TRUE) # counts of marked individuals unidentified
knownMarks <- matrix(covars[, "knownMarks"],
nrow = nGroups,
ncol = nocc,
byrow = TRUE) # counts of known marked individuals
effort <- matrix(covars[, "effort"],
nrow = nGroups,
ncol = nocc,
byrow = TRUE) # effort of observation in each colony at each timeNow you have all the basic data to run your models. But first, you need the models.
To define your models, you have to think of which variables will influence your parameters. For example: different effort in each occasion will influence the probability of observing individuals (α), so α could vary in function of effort or even, both effort and colony (‘colonies’) depending if you think that it is harder or easier to observe the individuals in different places. We can also come up with extra categories related to disturbance level (i.e. capim-açu is less disturbed than boldro), for example. Here I ran only one model so far, but you should later work on more models and we can develop a comparison test. :)
Basically you have to give a formula to each variable. for example: the basic formula formula = ~1 means that the parameter that this formula is assigned to is constant. We can also make this parameter vary with effort (formula = ~ effort) or even to effort and colonies as suggested before (formula = ~ effort + colonies). Within a model formula you can also set parameters to fixed values. For example, we know that no mabuyas died within a week of work, so we can set survival (Phi) to 1. We do that by using the constant formula and adding the argument fixed = 1 (formula = ~ 1, fixed = 1). It is also important that all models are in a list(). All parameters are descibed below in section 6. Make a data design.
alpha.effort <- list(formula = ~ effort)
Phi <- list(formula = ~ 1, fixed = 1)
U <- list(formula = ~ colonies)
constant <- list(formula = ~ 1) # Here I create a generic constant model to be used to terms that I need
zero <- list(formula = ~ 1, fixed = 0) # Here I create a generic zero model to be used to terms that I needTo have the data ready with the covariates you want to consider (i.e. effort), you need to go through 2 steps before running your model:
: use the function process.data()
: use the function make.design.data()
mabuya.process <- process.data(ch,
model = "PoissonMR",
time.intervals = rep(1, times = (nocc-1)),
groups = "colonies",
counts = list("Unmarked Seen" = unmarkedSeen,
"Marked Unidentified" = markedUnidentified,
"Known Marks" = knownMarks))
mabuya.process## $data
## ch colonies NA. group
## 1 010101010101 Americano NA 1
## 2 +0+001+0+0+0 Americano NA 1
## 3 +0+0+0+0+0+0 Americano NA 1
## 4 +0+0+0+0+0+0 Americano NA 1
## 5 +0+0+0+0+0+0 Americano NA 1
## 6 01+001+0+0+0 Americano NA 1
## 7 01010101+001 Americano NA 1
## 8 +0+0+001+0+0 Americano NA 1
## 9 +0+0+0+0+0+0 Americano NA 1
## 10 01+0+0+001+0 Americano NA 1
## 11 +001+0+0+0+0 Americano NA 1
## 12 +0+0+0+0+0+0 Americano NA 1
## 13 +0+0+0+0+0+0 Americano NA 1
## 14 +0+0+0+0+001 Americano NA 1
## 15 +0+0+0+0+0+0 Americano NA 1
## 16 +0+001+0+0+0 Americano NA 1
## 17 01+001+0+0+0 CapimAcu NA 2
## 18 01+0+0+0+001 CapimAcu NA 2
## 19 +001+0+001+0 CapimAcu NA 2
## 20 +0+0+0+00101 CapimAcu NA 2
## 21 +001+0+0+0+0 CapimAcu NA 2
## 22 010101010101 ForteBoldro NA 3
## 23 +0+00101+0+0 ForteBoldro NA 3
## 24 +0+0+0+0+0+0 ForteBoldro NA 3
## 25 +0+0+0+0+0+0 ForteBoldro NA 3
## 26 +00101+0+001 ForteBoldro NA 3
## 27 +0+00101+0+0 ForteBoldro NA 3
## 28 +0+0+001+001 ForteBoldro NA 3
## 29 +0+001+0+0+0 ForteBoldro NA 3
## 30 010101+0+0+0 ForteBoldro NA 3
## 31 +0+0+0+0+0+0 ForteBoldro NA 3
## 32 +0+001+0+0+0 ForteBoldro NA 3
## 33 +0+0+0+001+0 ForteBoldro NA 3
## 34 +001+001+0+0 ForteBoldro NA 3
## 35 +0+001+0+0+0 ForteBoldro NA 3
## 36 +0+0+0+0+001 Leao NA 4
## 37 +0+0+0010101 Leao NA 4
## 38 +001+0+0+0+0 Leao NA 4
## 39 +0+0+0+0+001 Leao NA 4
## 40 +0+0+0+0+0+0 Leao NA 4
## 41 +001+0+0+0+0 PedreiraSueste NA 5
## 42 +0+0+0+00101 PedreiraSueste NA 5
## 43 010101+0+0+0 PedreiraSueste NA 5
## 44 +0+0+0+0+0+0 PedreiraSueste NA 5
## 45 +0+0+0+00101 PedreiraSueste NA 5
## 46 +0+0+0+0+0+0 PedreiraSueste NA 5
## 47 01+0+001+001 Piquinho NA 6
## 48 0101010101+0 Piquinho NA 6
## 49 +0+0+001+0+0 Piquinho NA 6
## 50 0101+0+00101 Piquinho NA 6
## 51 01+0+0+0+0+0 PraiaBoldro NA 7
## 52 +0+0+0+0+0+0 PraiaBoldro NA 7
## 53 +0+0+0+0+0+0 PraiaBoldro NA 7
## 54 +0+0+0+0+0+0 PraiaBoldro NA 7
## 55 +0+0+0+0+0+0 PraiaBoldro NA 7
## 56 +0+0010101+0 PraiaBoldro NA 7
## 57 0101+0+001+0 PraiaBoldro NA 7
## 58 +0+0+0+0+0+0 PraiaBoldro NA 7
## 59 01+0+0+0+0+0 PraiaBoldro NA 7
## 60 01+0+0+0+001 PraiaBoldro NA 7
## 61 +0+0+001+0+0 PraiaBoldro NA 7
## 62 01+0+0+0+0+0 PraiaBoldro NA 7
## 63 01+0+0+0+0+0 PraiaBoldro NA 7
## 64 +0+0+0+0+0+0 PraiaBoldro NA 7
## 65 0101010101+0 PraiaBoldro NA 7
## 66 +00101010101 PraiaBoldro NA 7
## 67 +0+0+0+0+001 PraiaBoldro NA 7
## 68 +0+0+001+0+0 PraiaBoldro NA 7
## 69 +0+0+0+0+0+0 PraiaBoldro NA 7
## 70 +0+0+0+0+0+0 PraiaBoldro NA 7
## 71 +001010101+0 TejuAcu NA 8
## 72 +0+0+0+0+0+0 TejuAcu NA 8
## 73 +00101+00101 TejuAcu NA 8
##
## $model
## [1] "PoissonMR"
##
## $mixtures
## [1] 1
##
## $freq
## coloniesAmericano coloniesCapimAcu coloniesForteBoldro coloniesLeao
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 1 0 0 0
## 5 1 0 0 0
## 6 1 0 0 0
## 7 1 0 0 0
## 8 1 0 0 0
## 9 1 0 0 0
## 10 1 0 0 0
## 11 1 0 0 0
## 12 1 0 0 0
## 13 1 0 0 0
## 14 1 0 0 0
## 15 1 0 0 0
## 16 1 0 0 0
## 17 0 1 0 0
## 18 0 1 0 0
## 19 0 1 0 0
## 20 0 1 0 0
## 21 0 1 0 0
## 22 0 0 1 0
## 23 0 0 1 0
## 24 0 0 1 0
## 25 0 0 1 0
## 26 0 0 1 0
## 27 0 0 1 0
## 28 0 0 1 0
## 29 0 0 1 0
## 30 0 0 1 0
## 31 0 0 1 0
## 32 0 0 1 0
## 33 0 0 1 0
## 34 0 0 1 0
## 35 0 0 1 0
## 36 0 0 0 1
## 37 0 0 0 1
## 38 0 0 0 1
## 39 0 0 0 1
## 40 0 0 0 1
## 41 0 0 0 0
## 42 0 0 0 0
## 43 0 0 0 0
## 44 0 0 0 0
## 45 0 0 0 0
## 46 0 0 0 0
## 47 0 0 0 0
## 48 0 0 0 0
## 49 0 0 0 0
## 50 0 0 0 0
## 51 0 0 0 0
## 52 0 0 0 0
## 53 0 0 0 0
## 54 0 0 0 0
## 55 0 0 0 0
## 56 0 0 0 0
## 57 0 0 0 0
## 58 0 0 0 0
## 59 0 0 0 0
## 60 0 0 0 0
## 61 0 0 0 0
## 62 0 0 0 0
## 63 0 0 0 0
## 64 0 0 0 0
## 65 0 0 0 0
## 66 0 0 0 0
## 67 0 0 0 0
## 68 0 0 0 0
## 69 0 0 0 0
## 70 0 0 0 0
## 71 0 0 0 0
## 72 0 0 0 0
## 73 0 0 0 0
## coloniesPedreiraSueste coloniesPiquinho coloniesPraiaBoldro
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## 7 0 0 0
## 8 0 0 0
## 9 0 0 0
## 10 0 0 0
## 11 0 0 0
## 12 0 0 0
## 13 0 0 0
## 14 0 0 0
## 15 0 0 0
## 16 0 0 0
## 17 0 0 0
## 18 0 0 0
## 19 0 0 0
## 20 0 0 0
## 21 0 0 0
## 22 0 0 0
## 23 0 0 0
## 24 0 0 0
## 25 0 0 0
## 26 0 0 0
## 27 0 0 0
## 28 0 0 0
## 29 0 0 0
## 30 0 0 0
## 31 0 0 0
## 32 0 0 0
## 33 0 0 0
## 34 0 0 0
## 35 0 0 0
## 36 0 0 0
## 37 0 0 0
## 38 0 0 0
## 39 0 0 0
## 40 0 0 0
## 41 1 0 0
## 42 1 0 0
## 43 1 0 0
## 44 1 0 0
## 45 1 0 0
## 46 1 0 0
## 47 0 1 0
## 48 0 1 0
## 49 0 1 0
## 50 0 1 0
## 51 0 0 1
## 52 0 0 1
## 53 0 0 1
## 54 0 0 1
## 55 0 0 1
## 56 0 0 1
## 57 0 0 1
## 58 0 0 1
## 59 0 0 1
## 60 0 0 1
## 61 0 0 1
## 62 0 0 1
## 63 0 0 1
## 64 0 0 1
## 65 0 0 1
## 66 0 0 1
## 67 0 0 1
## 68 0 0 1
## 69 0 0 1
## 70 0 0 1
## 71 0 0 0
## 72 0 0 0
## 73 0 0 0
## coloniesTejuAcu
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
## 7 0
## 8 0
## 9 0
## 10 0
## 11 0
## 12 0
## 13 0
## 14 0
## 15 0
## 16 0
## 17 0
## 18 0
## 19 0
## 20 0
## 21 0
## 22 0
## 23 0
## 24 0
## 25 0
## 26 0
## 27 0
## 28 0
## 29 0
## 30 0
## 31 0
## 32 0
## 33 0
## 34 0
## 35 0
## 36 0
## 37 0
## 38 0
## 39 0
## 40 0
## 41 0
## 42 0
## 43 0
## 44 0
## 45 0
## 46 0
## 47 0
## 48 0
## 49 0
## 50 0
## 51 0
## 52 0
## 53 0
## 54 0
## 55 0
## 56 0
## 57 0
## 58 0
## 59 0
## 60 0
## 61 0
## 62 0
## 63 0
## 64 0
## 65 0
## 66 0
## 67 0
## 68 0
## 69 0
## 70 0
## 71 1
## 72 1
## 73 1
##
## $nocc
## [1] 6
##
## $nocc.secondary
## NULL
##
## $time.intervals
## [1] 1 1 1 1 1
##
## $begin.time
## [1] 1
##
## $age.unit
## [1] 1
##
## $initial.ages
## [1] 0 0 0 0 0 0 0 0
##
## $group.covariates
## colonies
## 1 Americano
## 2 CapimAcu
## 3 ForteBoldro
## 4 Leao
## 5 PedreiraSueste
## 6 Piquinho
## 7 PraiaBoldro
## 8 TejuAcu
##
## $nstrata
## [1] 1
##
## $strata.labels
## NULL
##
## $counts
## $counts$`Unmarked Seen`
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 17 9 17 8 12 13
## [2,] 2 3 3 1 1 0
## [3,] 8 9 17 13 13 12
## [4,] 1 3 2 0 0 3
## [5,] 2 2 2 3 2 3
## [6,] 4 5 7 4 6 4
## [7,] 22 35 32 30 28 23
## [8,] 1 3 3 3 1 3
##
## $counts$`Marked Unidentified`
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0
##
## $counts$`Known Marks`
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 16 16 16 16 16 16
## [2,] 5 5 5 5 5 5
## [3,] 14 14 14 14 14 14
## [4,] 5 5 5 5 5 5
## [5,] 6 6 6 6 6 6
## [6,] 4 4 4 4 4 4
## [7,] 20 20 20 20 20 20
## [8,] 3 3 3 3 3 3
##
##
## $reverse
## [1] FALSE
##
## $areas
## NULLNow you can see a list that summarizes your data. You can see the variables by typing names(mabuya.process):
## [1] "data" "model" "mixtures"
## [4] "freq" "nocc" "nocc.secondary"
## [7] "time.intervals" "begin.time" "age.unit"
## [10] "initial.ages" "group.covariates" "nstrata"
## [13] "strata.labels" "counts" "reverse"
## [16] "areas"And you can access these elements by typing mabuya.process$ELEMENT_NAME.
Now you need to make a data design that will identify the parameters in your data that are constant, and the variable ones for each model term, as well as set your covariates where they should be placed in:
mabuya.ddl <- make.design.data(mabuya.process)Using the names(mabuya.ddl) you can check now your parameters.
The , has have four main parameters:
(U) = number of unmarked individuals in the population during primary interval j. (alpha) = intercept (on log scale) for mean resighting rate during primary interval j. If there is no individual heterogeneity (= 0), once back-transformed from the log scale the real parameter estimate can be interpreted as the mean resighting rate for the entire population. (Phi) = apparent survival between primary intervals j and j + 1. (sigma) = individual heterogeneity level (on the log scale) during primary interval j, i.e., the additional variance due to a random individual heterogeneity effect with mean zero.
Now we can check one of the parameters: alpha:
mabuya.ddl$alphaNow you can include the covariates (such as effort and/or disturbance level) with the following code:
mabuya.ddl$alpha$effort <- as.numeric(apply(X = mabuya.ddl$alpha, 1, function(x){
col <- x[["group"]]
occ <- as.numeric(x[["time"]])
eff <- covars[covars$colonies==col&covars$occasion==occ, "effort"]
})
)And now we check again to see how this changed our data.frame:
mabuya.ddl$alphaHere, we added effort to alpha (which is the mean resighting rate during primary interval), as more effort increases the chances of observing an individual and the effort varied on occasions and colonies.
Now we are ready to run the model(s).
setwd(dir = file.path(getwd(), "outputs")) # Changing directory for outputs
mabuya.model <- mark(data = mabuya.process, # Here we use the processed data
ddl = mabuya.ddl, # Here we use the designed data
model.parameters = list( # Here are the parameters's models set on step 4.
alpha = alpha.effort,
Phi = Phi,
U = U,
sigma = zero,
GammaDoublePrime = constant,
GammaPrime = constant
)
)##
## Note: only 10 parameters counted of 12 specified parameters## AICc and parameter count have been adjusted upward##
## Output summary for PoissonMR model
## Name : alpha(~effort)sigma(~1)U(~colonies)Phi(~1)Gamma''(~1)Gamma'(~1)
##
## Npar : 12 (unadjusted=10)
## -2lnL: 708.574
## AICc : 733.2336 (unadjusted=729.03716)
##
## Beta
## estimate se lcl
## alpha:(Intercept) -1.3990600 0.1409582 -1.6753381
## alpha:effort 0.0019278 0.0045554 -0.0070009
## U:(Intercept) 3.8871439 0.1473813 3.5982766
## U:coloniesCapimAcu -2.0136670 0.3241977 -2.6490945
## U:coloniesForteBoldro -0.0649545 0.1619430 -0.3823628
## U:coloniesLeao -1.9003323 0.3436288 -2.5738448
## U:coloniesPedreiraSueste -1.9110580 0.2892956 -2.4780774
## U:coloniesPiquinho -0.9915231 0.2169601 -1.4167649
## U:coloniesPraiaBoldro 0.7962687 0.1369308 0.5278843
## U:coloniesTejuAcu -1.8249691 0.2900483 -2.3934638
## GammaDoublePrime:(Intercept) -22.6346940 4292.5022000 -8435.9392000
## GammaPrime:(Intercept) -0.1000000 1130.1295000 -2215.1539000
## ucl
## alpha:(Intercept) -1.1227819
## alpha:effort 0.0108564
## U:(Intercept) 4.1760112
## U:coloniesCapimAcu -1.3782395
## U:coloniesForteBoldro 0.2524538
## U:coloniesLeao -1.2268197
## U:coloniesPedreiraSueste -1.3440387
## U:coloniesPiquinho -0.5662813
## U:coloniesPraiaBoldro 1.0646530
## U:coloniesTejuAcu -1.2564744
## GammaDoublePrime:(Intercept) 8390.6699000
## GammaPrime:(Intercept) 2214.9539000
##
##
## Real Parameter alpha
## 1 2 3 4
## Group:coloniesAmericano 0.2718047 0.2540705 0.2590159 0.2565313
## Group:coloniesCapimAcu 0.2570263 0.2521189 0.2540705 0.2540705
## Group:coloniesForteBoldro 0.2640576 0.2540705 0.2590159 0.2590159
## Group:coloniesLeao 0.2555441 0.2565313 0.2530928 0.2550520
## Group:coloniesPedreiraSueste 0.2565313 0.2550520 0.2570263 0.2555441
## Group:coloniesPiquinho 0.2555441 0.2526054 0.2600165 0.2540705
## Group:coloniesPraiaBoldro 0.2555441 0.2625349 0.2666151 0.2605182
## Group:coloniesTejuAcu 0.2540705 0.2506650 0.2506650 0.2516333
## 5 6
## Group:coloniesAmericano 0.2540705 0.2540705
## Group:coloniesCapimAcu 0.2540705 0.2560372
## Group:coloniesForteBoldro 0.2625349 0.2590159
## Group:coloniesLeao 0.2530928 0.2540705
## Group:coloniesPedreiraSueste 0.2545608 0.2535812
## Group:coloniesPiquinho 0.2540705 0.2516333
## Group:coloniesPraiaBoldro 0.2526054 0.2666151
## Group:coloniesTejuAcu 0.2516333 0.2516333
##
##
## Real Parameter sigma
## 1 2 3 4 5 6
## Group:coloniesAmericano NA NA NA NA NA NA
## Group:coloniesCapimAcu NA NA NA NA NA NA
## Group:coloniesForteBoldro NA NA NA NA NA NA
## Group:coloniesLeao NA NA NA NA NA NA
## Group:coloniesPedreiraSueste NA NA NA NA NA NA
## Group:coloniesPiquinho NA NA NA NA NA NA
## Group:coloniesPraiaBoldro NA NA NA NA NA NA
## Group:coloniesTejuAcu NA NA NA NA NA NA
##
##
## Real Parameter U
## 1 2 3 4
## Group:coloniesAmericano 48.771390 48.771390 48.771390 48.771390
## Group:coloniesCapimAcu 6.510895 6.510895 6.510895 6.510895
## Group:coloniesForteBoldro 45.704164 45.704164 45.704164 45.704164
## Group:coloniesLeao 7.292246 7.292246 7.292246 7.292246
## Group:coloniesPedreiraSueste 7.214449 7.214449 7.214449 7.214449
## Group:coloniesPiquinho 18.094731 18.094731 18.094731 18.094731
## Group:coloniesPraiaBoldro 108.138470 108.138470 108.138470 108.138470
## Group:coloniesTejuAcu 7.863052 7.863052 7.863052 7.863052
## 5 6
## Group:coloniesAmericano 48.771390 48.771390
## Group:coloniesCapimAcu 6.510895 6.510895
## Group:coloniesForteBoldro 45.704164 45.704164
## Group:coloniesLeao 7.292246 7.292246
## Group:coloniesPedreiraSueste 7.214449 7.214449
## Group:coloniesPiquinho 18.094731 18.094731
## Group:coloniesPraiaBoldro 108.138470 108.138470
## Group:coloniesTejuAcu 7.863052 7.863052
##
##
## Real Parameter Phi
## Group:coloniesAmericano
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesCapimAcu
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesForteBoldro
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesLeao
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesPedreiraSueste
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesPiquinho
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesPraiaBoldro
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
## Group:coloniesTejuAcu
## 1 2 3 4 5
## 1
## 2
## 3
## 4
## 5
##
##
## Real Parameter GammaDoublePrime
## Group:coloniesAmericano
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesCapimAcu
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesForteBoldro
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesLeao
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesPedreiraSueste
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesPiquinho
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesPraiaBoldro
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
## Group:coloniesTejuAcu
## 1 2 3 4 5
## 1 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 2 1.47869e-10 1.47869e-10 1.47869e-10 1.47869e-10
## 3 1.47869e-10 1.47869e-10 1.47869e-10
## 4 1.47869e-10 1.47869e-10
## 5 1.47869e-10
##
##
## Real Parameter GammaPrime
## Group:coloniesAmericano
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesCapimAcu
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesForteBoldro
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesLeao
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesPedreiraSueste
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesPiquinho
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesPraiaBoldro
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208
##
## Group:coloniesTejuAcu
## 2 3 4 5
## 1 0.4750208 0.4750208 0.4750208 0.4750208
## 2 0.4750208 0.4750208 0.4750208
## 3 0.4750208 0.4750208
## 4 0.4750208setwd(dir = directory) # Changing directory back to originalNow, it is important to notice that the U the model outputs is NOT the population size. The last, as well as other parameters can, however be derived:
Derived Parameters
= overall mean resighting rate for primary occasion j. This is a parameter derived as a function of , ^2[j] and and E[j]. Note that when = 0 and E[j] = 0, then the real parameter estimate for is identical to the derived parameter estimate for .
N[j] = U[j] + n[j], or the total population size during primary occasion j. This is a derived parameter, because MARK actually estimates U[j] in the model. n[j] is the marked population.
RMark automatically calculates the population size for each group (and occasion when we do a long term study with longer interval time).
mabuya.model$results$derived$`N Population Size`You will notice, however, that it doesn’t output the location’s names. We can, however retrieve it from the original data, and clean it up (repetitions):
populationSize <- cbind(covars$colonies, mabuya.model$results$derived$`N Population Size`) %>%
unique()
populationSize <- populationSize[populationSize$se > 0,]Now we can also calculate the densities among the areas (and at some point compare them statistically):
populationSize$areaSize <- data.frame(areaSize = c(36.7, 44.92, 7.7, 13.5, 15, 52, 48, 30))
populationSize$density <- data.frame(density = populationSize$estimate/populationSize$areaSize)
names(populationSize) <- c("locations", "estimate", "standard.error", "lowerCI", "upperCI", "areaSize", "density")And here we see the final result of the population size: