Vamos usar o pacote RMark, que é um pacote do R para usar o programa MARK. Siga as instruções do site do RMark para instalar o pacote.
Com o RMark instalado, abra o R e carregue-o:
library(RMark)Usaremos dados de registros fotográficos de indivíduos do boto cinza (Sotalia guianensis) em 11 ocasiões. Aqui há mais informações sobre este caso de estudo.
Os dados estão no formato nativo do MARK (.inp). Use os comandos abaixo para importá-lo para o R:
## Link dos dados na página da disciplina
url <- "http://ecologia.ib.usp.br/bie5703/lib/exe/fetch.php?media=roteiros:botos_2002.inp"
## Importa arquivo inp
boto2002 <- convert.inp(url)O primeiro passo é usar a função process.data para criar um objeto com as informações que o Mark usa para ajustar o modelo. Uma delas é o tipo de modelo, que é indicado no argumento model.
Para o modelo de populações fechada sem heterogeneidade e de verossimilhança condicionada este argumento é model="Closed" 1:
boto <- process.data(data=boto2002, model="Closed")E para o modelo com heterogeneidade o argumento é model="FullHet"
botoH <- process.data(data=boto2002, model="FullHet")Para ajustar os modelos, crie listas que especificam a fórmula de cada termo. No modelo Closed os nomes parâmetros que podem variar são p (\(p\) ,probabilidade da primeira captura) e c (\(c\), probabilidade de recaptura) 2. O objeto criado na seção acima tem uma covariável de tempo chamada time, que então pode ser usado nas fórmulas para expressar diferenças entre ocasiões 3:
## Fórmulas estatísticas para cada parâmetro do modelo sem heterogeneidade
## p e c constantes mas diferentes
t.dot <- list(formula=~1)
## p=c contantes (use o argumento share=TRUE)
t.dotshared=list(formula=~1,share=TRUE)
## Parametros dependem do tempo
t.time <- list(formula=~time)
## Parametro p=c dependem do tempo
t.timeshared <- list(formula=~time, share=TRUE)E usamos a função mark para fazer os ajuste 4:
boto.M0 <- mark(boto, model.parameters=list(p=t.dotshared),
model.name="M0", adjust=TRUE, delete=TRUE)
boto.Mb <- mark(boto, model.parameters=list(c=t.dot, p=t.dot),
model.name="Mb", adjust=TRUE, delete=TRUE)##
## Note: only 2 parameters counted of 3 specified parameters
##
## AICc and parameter count have been adjusted upwardboto.Mt <- mark(boto, model.parameters=list(p=t.timeshared), model.name="Mt",
adjust=TRUE, delete=TRUE)
boto.Mtb <- mark(boto, model.parameters=list(c=t.time, p=t.time),
model.name="Mtb", adjust=TRUE, delete=TRUE)##
## Note: only 18 parameters counted of 22 specified parameters
##
## AICc and parameter count have been adjusted upwardSe omitimos a função de um parâmetro ela será constante. Portanto para todos os modelos acima a expressão para o parâmetro \(f_0\) é formula=~1.
Para os modelos com heterogeneidade acrescente o termo mixture nas fórmulas do parâmetro \(p\). O default é uma mistura de duas subpopulações, o que representa que uma proporção \(\pi\) dos indivíduos tem uma probabilidade de captura/recaptura e o restante (\(1-\pi\)) tenha outra.
No modelo com efeito de ocasião, use uma fórmula com interação. Com isso as probabilidades de captura e recaptura de cada subpopulação poderão ser diferentes a cada ocasião.
## Fórmulas estatísticas para cada parâmetro do modelo com heterogeneidade
## p com heterogeneidade
t.mix <- list(formula=~mixture)
## p=c com heterogeneidade (use o argumento share=TRUE)
t.mixshared=list(formula=~mixture,share=TRUE)
## Parametros diferem entre ocasiões
t.timemixshared <- list(formula=~time*mixture, share=TRUE)
t.timemix <- list(formula=~time*mixture)E ajuste os modelos
boto.Mh <- mark(botoH, model.parameters=list(p=t.mixshared),
model.name="Mh", adjust=TRUE, delete=TRUE)
boto.Mbh <- mark(botoH, model.parameters=list(c=t.mix, p=t.mix),
model.name="Mbh", adjust=TRUE, delete=TRUE)
boto.Mth <- mark(botoH, model.parameters=list(p=t.timemixshared),
model.name="Mth", adjust=TRUE, delete=TRUE)##
## Note: only 19 parameters counted of 24 specified parameters
##
## AICc and parameter count have been adjusted upwardboto.Mtbh <- mark(botoH, model.parameters=list(c=t.timemix, p=t.timemix),
model.name="Mtbh", adjust=TRUE, delete=TRUE)##
## Note: only 25 parameters counted of 44 specified parameters
##
## AICc and parameter count have been adjusted upwardNote que em todos os modelos acima os parâmetros \(f_0\) e \(\pi\) são constantes, pois omitimos suas fórmulas.
A função abaixo retorna a tabela de seleção de modelos:
collect.models(lx=c("boto.M0", "boto.Mb", "boto.Mt", "boto.Mtb",
"boto.Mh", "boto.Mbh", "boto.Mth", "boto.Mtbh"))## Warning in model.table(x, type, pf = 2, adjust = adjust): Model list contains models of differing types## model npar AICc DeltaAICc weight Deviance
## 7 Mth 24 252.9300 0.00000 9.997964e-01 148.1513
## 3 Mt 12 270.1989 17.26893 1.778327e-04 191.7698
## 6 Mbh 6 275.2837 22.35377 1.399112e-05 209.4365
## 5 Mh 4 275.6888 22.75878 1.142628e-05 213.9520
## 4 Mtb 22 282.6256 29.69561 3.561171e-07 182.3529
## 1 M0 2 289.7349 36.80490 1.018190e-08 232.0679
## 2 Mb 3 290.1108 37.18087 8.436994e-09 230.4140
## 8 Mtbh 44 291.2933 38.36331 4.671170e-09 138.7168A função coef retorna os coeficientes na escala de ligação (logito). Para as estimativas na escala de probabilidades use a função get.real:
coef(boto.Mth, data=boto2002)## estimate se lcl ucl
## pi:(Intercept) -1.9020436 0.4820132 -2.846789 -0.9572976
## p:(Intercept) 18.9949060 130.1367200 -236.073060 274.0628800
## p:time2 -18.5895140 130.1395100 -273.662950 236.4839200
## p:time3 -18.5895140 130.1395100 -273.662950 236.4839200
## p:time4 -17.6086330 130.1400300 -272.683100 237.4658400
## p:time5 13.4694670 505.2418300 -976.804540 1003.7435000
## p:time6 -19.4004230 130.1394700 -274.473800 235.6729500
## p:time7 -18.5894250 130.1390000 -273.661870 236.4830200
## p:time8 -19.4003390 130.1398300 -274.474420 235.6737400
## p:time9 14.7536400 0.1681218 14.424121 15.0831590
## p:time10 -19.4003350 130.1402300 -274.475200 235.6745300
## p:time11 -19.4003650 130.1402100 -274.475180 235.6744500
## p:mixture2 -20.3260450 130.1368100 -275.394190 234.7421000
## p:time2:mixture2 18.9193500 130.1401400 -236.155320 273.9940200
## p:time3:mixture2 18.9193520 130.1401400 -236.155320 273.9940300
## p:time4:mixture2 18.6085240 130.1405100 -236.466880 273.6839300
## p:time5:mixture2 -12.7213890 505.2418400 -1002.995400 977.5526300
## p:time6:mixture2 19.7302310 130.1401000 -235.344380 274.8048400
## p:time7:mixture2 19.3375060 130.1395100 -235.735940 274.4109500
## p:time8:mixture2 -26.3891640 0.0522470 -26.491568 -26.2867600
## p:time9:mixture2 -45.6317130 0.0874371 -45.803090 -45.4603370
## p:time10:mixture2 18.4124160 130.1418100 -236.665540 273.4903700
## p:time11:mixture2 18.7334550 130.1414500 -236.343790 273.8107000
## f0:(Intercept) 0.4039731 1.1145751 -1.780594 2.5885403## Na escala de probabilidades
get.real(boto.Mth, parameter="p")## [[1]]
## 1 2 3 4 5 6
## mixture:1 1.000000 0.5999825 0.5999825 0.7999967 1.0000000 0.3999876
## mixture:2 0.208971 0.2686853 0.2686856 0.4179369 0.3582285 0.2686798
## 7 8 9 10 11
## mixture:1 0.6000038 4.000078e-01 1.000000e+00 0.4000088 0.4000016
## mixture:2 0.3582291 3.433721e-21 1.027346e-14 0.0895568 0.1194079get.real(boto.Mth, parameter="pi")## [[1]]
##
## mixture:1 0.1298774get.real(boto.Mth, parameter="f0") ##N de indivíduos não registrados## [[1]]
## 1
## 1.497764A estimativa de interesse é o tamanho da população, que é obtido somando-se a \(f_0\) ao total de indivíduos registrados. Usamos a função get.real com argumento se=TRUE para obter os intervalos de confiança5:
(boto.f0.ic <- as.numeric(get.real(boto.Mth, parameter="f0", se=TRUE)[,5:6]))## [1] 0.2573562 8.7166958O número de indivíduos registrados é a soma das frequências no objeto processado
(boto.Nobs <- sum(boto$freq))## [1] 37E finalmente temos o intervalo de confiança do tamanho populacional
(boto.Nobs + boto.f0.ic)## [1] 37.25736 45.71670A estimativa do tamanho populacional parece bastante precisa, mas com o pacote Rcapture os intervalos são mais conservadores. Confira isto executando o roteiro do Rcapture
A lista de modelos implementados no RMark está no diretório onde o R instalou o pacote. Você pode consultá-lo lá ou no repositório de desenvolvimento do RMark: (https://github.com/jlaake/RMark/blob/master/RMark/inst/MarkModels.pdf).↩
O parâmetro \(f_0\) (número de indivíduos não registrados em nenhuma ocasião) é constante por definição, já que a população é fechada.↩
Para entender completamente isso estude o comando make.design.data e o objeto que ele cria, que é a uma lista de matrizes de delineamento do modelo. Se você entender este objeto saberá quais covariáveis estão disponíveis e como manipulá-las. Veja também o apêndice sobre o RMark no guia online do MARK.↩
Use sempre os argumentos delete=TRUE para remover os arquivos temporários do Mark que o ajuste cria e adjust=TRUE para que o número de parâmetros seja verificado e ajustado, para o cálculo do AIC.↩
Para uma lista com todos os coeficientes e seus intervalos use a função summary com a opção se=TRUE.↩