Как бороться с временной корреляцией / трендом в миллионных долях в минуту

Добрый день,

Я работал с Baddeley et al. 2015, чтобы подогнать модель точечного процесса к нескольким точечным образцам с помощью < em> mppm {spatstat}. Мои точечные шаблоны - это годовые данные о количестве крупных травоядных (то есть точечные участки (x, y) самцов / самок * 5 лет) на охраняемой территории (owin). У меня есть несколько пространственных ковариат, например. расстояние до рек (rivD) и продуктивность растительности (NDVI).

Первоначально я подобрал модель, в которой реакция травоядных была функцией rivD + NDVI, и позволил коэффициентам варьироваться в зависимости от пола (см. Mppm1 в воспроизводимом примере ниже). Однако мои годовые точечные диаграммы не являются независимыми между годами в том смысле, что существует тенденция к увеличению во времени (т.е. в год 1 экспоненциально больше животных по сравнению с годом 5). Поэтому я добавил год как случайный эффект, думая, что если я позволю перехвату изменяться за год, я смогу это учесть (см. Mppm2).

Теперь мне интересно, правильный ли это способ сделать это? Если бы я использовал GAMM gamm {mgcv}, я бы добавил структуру временной корреляции, например correlation = corAR1(form=~year) но не думаете, что это возможно в mppm (см. Mppm3)?

Я был бы очень признателен за любые идеи о том, как работать с этой структурой временной корреляции в реплицированном точечном шаблоне с mppm {spatstat}.

Большое тебе спасибо

Сандра

# R version 3.3.1 (64-bit)
library(spatstat) # spatstat version 1.45-2.008

#### Simulate point patterns
# multitype Neyman-Scott process (each cluster is a multitype process)
nclust2 = function(x0, y0, radius, n, types=factor(c("male", "female"))) {
  X = runifdisc(n, radius, centre=c(x0, y0))
  M = sample(types, n, replace=TRUE)
  marks(X) = M
  return(X)
}

year1 = rNeymanScott(5,0.1,nclust2, radius=0.1, n=5)
# plot(year1)
#-------------------
year2 = rNeymanScott(10,0.1,nclust2, radius=0.1, n=5)
# plot(year2)
#-------------------
year2 = rNeymanScott(15,0.1,nclust2, radius=0.1, n=10)
# plot(year2)
#-------------------
year3 = rNeymanScott(20,0.1,nclust2, radius=0.1, n=10)
# plot(year3)
#-------------------
year4 = rNeymanScott(25,0.1,nclust2, radius=0.1, n=15)
# plot(year4)
#-------------------
year5 = rNeymanScott(30,0.1,nclust2, radius=0.1, n=15)
# plot(year5)

#### Simulate distance to rivers
line <- psp(runif(10), runif(10), runif(10), runif(10), window=owin())
# plot(line)
# plot(year1, add=TRUE)

#------------------------ UPDATE ------------------------#
#### Create hyperframe
#---> NDVI simulated with distmap to point patterns (not ideal but just to test)
hyp.years = hyperframe(year=factor(2010:2014),
                       ppp=list(year1,year2,year3,year4,year5),
                       NDVI=list(distmap(year5),distmap(year1),distmap(year2),distmap(year3),distmap(year4)),
                       rivD=distmap(line),
                       stringsAsFactors=TRUE)
hyp.years$numYear = with(hyp.years,as.numeric(year)-1)
hyp.years

#### Run mppm models
# mppm1 = mppm(ppp~(NDVI+rivD)/marks,data=hyp.years); summary(mppm1)
#..........................
# mppm2 = mppm(ppp~(NDVI+rivD)/marks,random = ~1|year,data=hyp.years); summary(mppm2)
#..........................
# correlation = corAR1(form=~year)
# mppm3 = mppm(ppp~(NDVI+rivD)/marks,correlation = corAR1(form=~year),use.gam = TRUE,data=hyp.years); summary(mppm3)

###---> Run mppm model with annual trend and random variation in growth
mppmCorr = mppm(ppp~(NDVI+rivD+numYear)/marks,random = ~1|year,data=hyp.years)
summary(mppm1)

person Sands    schedule 26.06.2016    source источник


Ответы (1)


Если наблюдается тренд в размере популяции с течением времени, то, возможно, имеет смысл включить эту тенденцию в систематическую часть модели. Я бы посоветовал вам добавить новую числовую переменную NumYear во фрейм данных (например, указав количество лет, прошедших с 2010 года). Затем попробуйте добавить простые термины тренда, такие как +NumYear, в формулу модели (это будет соответствовать экспоненциальному росту популяции, который вы наблюдали). Вы можете оставить 1|year член случайного эффекта, который затем учтет случайные вариации в размере популяции около долгосрочная тенденция роста.

Нет необходимости разделять шаблоны данных для каждого года на отдельные мужские и женские модели. Переменная marks в формуле модели может использоваться для указания любой модели, которая зависит от пола.

Я почти уверен, что mppm с use.gam=TRUE не распознает аргумент correlation, и это, вероятно, просто игнорируется. (Это зависит от того, что происходит внутри gam).

person Adrian Baddeley    schedule 27.06.2016
comment
Спасибо, профессор Баддели за вашу помощь. Я обновил свой сценарий выше, как вы предложили. У меня разделились самец и самка, потому что я не мог получить оценки для работы в формуле. Я пробовал использовать обе метки и «метки», но продолжаю получать следующую ошибку: неверная формула модели в ExtractVars. Какие-либо предложения? - person Sands; 27.06.2016
comment
Спасибо, что заметили это. Это должно быть ошибка. Вы согласны с @ adrian-baddeley? Я внес небольшую поправку в код mppm, который, надеюсь, решит проблему с отметками в формуле. Вы должны использовать marks без кавычек. Вы можете протестировать эту новую версию, установив временную ветвь spatstat под названием mark_mppm (требуется devtools пакет): devtools::install_github("spatstat/spatstat", ref = "marked_mppm") - person Ege Rubak; 27.06.2016
comment
Превосходно! Большое тебе спасибо. Я соответствующим образом обновил приведенный выше сценарий. - person Sands; 27.06.2016
comment
Черт возьми. Да, это ошибка. Исправление Эге правильное. [Исправление ошибки будет включено в основную ветку разрабатываемой версии spatstat.] - person Adrian Baddeley; 28.06.2016