Байесовская модель смешанных эффектов для повторных измерений между годами в WINBUGS - PullRequest
Новая
       45

Байесовская модель смешанных эффектов для повторных измерений между годами в WINBUGS

0 голосов
/ 09 апреля 2019

Я заинтересован в подборе байесовской иерархической модели в WinBUGS, которая имеет дело с повторными измерениями для каждого местоположения в течение ряда лет.

Я смоделировал некоторые данные, чтобы показать, что я пытаюсь подогнать.Всего два года, и каждый год посещается несколько мест.В 1-м году было посещено 5 мест, а во 2-м году также было посещено 2 новых места, в которых было посещено 7 мест.Каждое местоположение имеет соответствующие значения еженедельного подсчета (недели = 15).

Еженедельные измерения определяются как разность функций между неделями.Где функция определяется как: 

f(week) = alpha\*exp(-exp(beta - gamma\*(week - delta))) - alpha\*exp(-exp(beta - gamma\*((week-1) - delta)))

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

library(R2WinBUGS)
set.seed(902010)

#define Modified Gompertz
Gompertz_fn <- function(alpha, beta, gamma, delta, t){
  round(alpha*exp(-exp(beta-gamma*(t-delta))))
}

#coefficients for GC, with alpha and beta different between years
alpha <- list(alpha1 = round(rnorm(5,100,1)), alpha2=round(rnorm(7, 75, 1))) 
beta <-  list(beta1 = (rnorm(5,1.5,0.05)), beta2=(rnorm(7,1, 0.05)))
gamma <- 0.5
delta<- 2.5
weeks <- 1:15

#matrix to store observed values
y1 <- matrix(NA,nrow=15,ncol=5) #year one, five areas for year one, 15 measurements
y2 <- matrix(NA,nrow=15,ncol=7) #year two, seven areas for year two, 15 measurements

#obtain observed counts for each area for year 1 and 2
for(i in 1:5){
y1[,i] <- Gompertz_fn(alpha=alpha[[1]][i], beta=beta[[1]][i], gamma=gamma, delta=delta, t=weeks) }

for(i in 1:7){
y2[,i] <- Gompertz_fn(alpha=alpha[[2]][i], beta=beta[[2]][i], gamma=gamma, delta=delta, t=weeks) }


#save in data frame
df <- data.frame(obs= c(y1,y2), area=c(rep(1:5,each=15), rep(1:7,each=15)), 
                 week=rep(1:15,12), year = c(rep(1, 5*15), rep(2,7*15)))

#Winbugs model
##########
sink("model1.txt")
cat("
    model{

    for(i in 1:n.obs){ #loop through observations

    y[i] ~ dpois(mu_area[i]) 

    log(mu_area[i]) <- alpha[area[i]]*exp(-exp(beta[area[i]] - gamma[year[i]]*(week[i] - delta[year[i]]))) - 

    alpha[area[i]]*exp(-exp(beta[area[i]] - gamma[year[i]]*((week[i]-1) - delta[year[i]])))
    }


    for(i in 1:n.areas){ #loop through areas
    alpha[area.2[i]] ~ dnorm(mu_alpha[years_areas[i]],tau_alpha[years_areas[i]])
    beta[area.2[i]] ~ dnorm(mu_beta[years_areas[i]],tau_beta[years_areas[i]])
    }

    for(i in 1:n.years){ #loop through years
    mu_alpha[i] ~dnorm(0,0.001)
    mu_beta[i] ~dnorm(0,0.001)

    tau_alpha[i]~dgamma(0.1,1)
    tau_beta[i] ~ dgamma(0.1,1)

    gamma[i] ~ dnorm(0,0.1)
    delta[i] ~ dnorm(0,0.1)

    }



    }
    ",fill=TRUE)
sink()
##############

##################################################################################
#winbugs data
df1 <- list(y=df$obs,year= df$year,
            week=df$week,area=df$area,
            n.obs=180,n.areas=12,n.years=2 ,
            years_areas=c(rep(1,5), rep(2,7)),
            area.2 = 1:12)

#params to save
params <- c("alpha","beta","gamma","delta",
            "mu_alpha","mu_beta","tau_alpha","tau_beta")

#initial values
inits <- list(alpha=c(rep(100,5),rep(75,7)), 
              beta = c(rep(3,5), rep(2,7)), 
              gamma=rep(0.5,2), delta=rep(0,2),
              mu_alpha=c(100,50), mu_beta=c(3,2),
              tau_alpha=rep(1,2), tau_beta=rep(1,2))

#run model
b2 <- bugs(data = df1, inits = list(inits),
           parameters.to.save = params,
           n.chains = 1, n.iter = 10000, 
           n.burnin =5000,
           n.thin = 10, DIC=T,
           codaPkg = F,
           model.file = "model1.txt", 
           bugs.directory = paste0(Sys.getenv(c("USERPROFILE")), "\\WinBUGS14"),debug=T)

Я ожидаю получить параметры альфа и бета для каждого местоположения и года.гамма, дельта, mu_alpha, tau_alpha, mu_beta и tau_beta для конкретного года.Однако мои результаты показывают, что альфа и бета зависят от области, но не различаются по годам.

...