fpemmodeling: fpemmodeling

write_jags_model <- function(old_dm = FALSE,
                             include_ss_data = FALSE,
                             nulldata,
                             is_in_union) {
  
  # model taken from https://github.com/FPcounts/ContraceptiveUse_MWRA_Update
  # deleting all global-only stuff
  # remove jumping over years w/o observations
  # kept C because easier
  # but code is coded for C= 1 (ie hardcoded subreg index)
  
  cat("
model{
  Y  ~ dnorm(greghack, 1)
  greghack ~ dnorm(1,1)
  # AR part
  for (c in 1:C){
    # eps.P.ct[c, 1] ~ dnorm(0, tau.P.st) # for P
    # eps.R.ct[c,1] ~ dnorm(0, tau.R.st) # for R
    # eps.Z.ct[c,1] ~ dnorm(0, tau.Z.st) # for Z
    # for (t in 2:nyears){
    #   eps.P.ct[c, t] ~ dnorm(rho.P*eps.P.ct[c, t-1], tau.P)
    #   eps.R.ct[c, t] ~ dnorm(rho.R*eps.R.ct[c, t-1], tau.R)
    #   eps.Z.ct[c, t] ~ dnorm(rho.Z*eps.Z.ct[c, t-1], tau.Z)
    # }
    eps.P.ct[c, t_star] ~ dnorm(0, tau.P.st) # for P
    eps.R.ct[c, t_star] ~ dnorm(0, tau.R.st) # for R
    eps.Z.ct[c, t_star] ~ dnorm(0, tau.Z.st) # for Z
    for (t in (t_star+1):nyears){
      eps.P.ct[c, t] ~ dnorm(rho.P*eps.P.ct[c, t-1], tau.P)
      eps.R.ct[c, t] ~ dnorm(rho.R*eps.R.ct[c, t-1], tau.R)
      eps.Z.ct[c, t] ~ dnorm(rho.Z*eps.Z.ct[c, t-1], tau.Z)
    }
    for (t in 2:t_star){
      eps.P.ct[c, t-1] ~ dnorm(rho.P*eps.P.ct[c, t], tau.P)
      eps.R.ct[c, t-1] ~ dnorm(rho.R*eps.R.ct[c, t], tau.R)
      eps.Z.ct[c, t-1] ~ dnorm(rho.Z*eps.Z.ct[c, t], tau.Z)
    }
  }
  tau.P.st <- tau.P*(1-pow(rho.P,2))
  tau.P <- pow(sigma.P, -2)
  tau.R.st <- tau.R*(1-pow(rho.R,2))
  tau.R <- pow(sigma.R, -2)
  tau.Z.st <- tau.Z*(1-pow(rho.Z,2))
  tau.Z <- pow(sigma.Z, -2)
  # logistic curves and model for Z
  for (c in 1:C){
    for (t in 1:nyears){
      Rmu.ct[c, t] <- Romega.c[c]*(t - RT.c[c])
      Rstar.ct[c, t] <- Rmax.c[c]/(1+exp(-Rmu.ct[c, t]))
      R.ct[c, t] <- 1/(1+exp(-( logit(Rstar.ct[c,t]) + eps.R.ct[c,t])))
      logitZstar.ct[c,t] <- (unmet.intercept.c[c]
                             + a.unmet
                             + b.unmet * (P.ct[c,t] - pmid.for.unmet)
                             + c.unmet * pow(P.ct[c,t] - pmid.for.unmet,2))
      Z.ct[c,t] <- 1/(1+exp(-(logitZstar.ct[c,t] + eps.Z.ct[c,t])))
      #neg.explogitZ.ct[c,t] = exp(-logitZ.ct[c,t])
    }
    for(t in 1:(t_star-1)){
      ls.ct[c,(t_star-t)] <- s.ct[c, (t_star-t)+1] - eps.P.ct[c, t_star-t] #logit
      ils.ct[c,(t_star-t)] <- 1/(1+exp(-ls.ct[c,(t_star-t)])) #inv.logit
      #Step function; test for x >/= 0
      I[c,(t_star-t)] <- step(ils.ct[c,(t_star-t)] - pmax.c[c])
      ###Get P.ct directly in the backward direction
      #Only need this bit if I=0 i.e., ils.ct<pmax.c
      zeta.ct[c,(t_star-t)] <- (1-I[c,(t_star-t)])*(logit(min((1-0.00001),ils.ct[c,(t_star-t)]/pmax.c[c]))-omega.c[c])
      P.ct[c,(t_star-t)]<-(1-I[c,(t_star-t)])*(pmax.c[c]*(1/(1+exp(-zeta.ct[c,(t_star-t)])))) + I[c,(t_star-t)]*ils.ct[c,(t_star-t)]
      ###Get logit(P.ct)
      s.ct[c,(t_star-t)] <- logit(P.ct[c,(t_star-t)])
    } # end back extrapolation
    for(t in (t_star+1):nyears){
      #Step function; test for x >/= 0
      I[c,t] <- step(P.ct[c,t-1] - pmax.c[c])
      #Only need this bit if I=0 i.e., P.ct<pmax.c
      zeta.ct[c,t] <- (1-I[c,t])*(logit(min((1-0.000001),P.ct[c,t-1]/pmax.c[c])) + omega.c[c])
      s.ct[c,t] <- logit(I[c,t]*(P.ct[c,t-1]) + (1-I[c,t])*pmax.c[c]*(1/(1+exp(-zeta.ct[c,t])))) + eps.P.ct[c,t-1]
      P.ct[c,t] <- 1/(1 + exp(-s.ct[c,t]))
    }
    ### add pmax_lower_bound here
    pmax.c[c] <- pmax_lower_bound + (1-pmax_lower_bound)/(1+exp(-logitpmax.c[c]))
    logitpmax.c[c] ~ dnorm(lp.world, 1/(pow(sigma.lpc, 2) + pow(sd_lp.world, 2)))
# lower bound for rmax is 0.5 for married AND unmarried
    Rmax.c[c] <- 0.5 + (1-0.5)/(1+exp(-logitRmax.c[c]))
    logitRmax.c[c] ~ dnorm(lr.world, 1/(pow(sigma.lrc, 2) + pow(sd_lr.world, 2)))
    logitomega.c[c] ~ dnorm(w.subreg, 1/(pow(sigma.wc, 2) + pow(sd_w.subreg, 2)))
    omega.c[c] <- 0.01 + (0.5-0.01)/(1+exp(-logitomega.c[c]))
    Romega.c[c] <- 0.01 + (0.5-0.01)/(1+exp(-logitRomega.c[c]))
    logitRomega.c[c] ~ dnorm(Rw.subreg, 1/(pow(sigma.Rwc, 2) + pow(sd_Rw.subreg, 2)))
    s.ct[c,t_star] <- setlevel.c[c]
    setlevel.c[c] ~ dnorm(mean_setlevel, 1/(var_setlevel + pow(sd_mean_setlevel, 2)))
    P.ct[c,t_star] <- 1/(1+exp(-s.ct[c,t_star]))
    RT.c[c] ~ dnorm(RT.subreg, 1/(pow(sigma.RTc, 2) + pow(sd_RT.subreg, 2)))
    unmet.intercept.c[c] ~ dnorm(unmet.subreg, 1/(pow(sigma.unmetc, 2) + pow(sd_unmet.subreg, 2)))
  } # end country loop
  # tau.lrc <- pow(sigma.lrc,-2)
  # tau.lpc <- pow(sigma.lpc,-2)
  # tau.wc <- pow(sigma.wc, -2)
 # tau.Sc <- pow(sigma.Sc,-2)
 # tau.higherSc <- pow(sigma.Sc,-2)
  # tau.Rwc <- pow(sigma.Rwc, -2)
  # tau.RTc <- pow(sigma.RTc, -2)
  # tau.unmetc <- pow(sigma.unmetc,-2)
  #------
  # to export
  for (c in 1:C){
    for (t in 1:nyears){
      mod.ct[c,t] <- P.ct[c,t]*R.ct[c,t]
      trad.ct[c,t] <- P.ct[c,t]*(1-R.ct[c,t])
      unmet.ct[c,t] <- (1-P.ct[c,t])*Z.ct[c,t]
    }
  }
  ",sep="",append=FALSE, file = "model.txt", fill = TRUE)
  
  # DMS
  # update name of observed props to avoid confusion
  if (!old_dm & !nulldata){ # simple
    ## note: does not work yet for partial/all missing
    cat("
    for (j in 1:n_mod){
      modern.j[j] = mod.ct[get_c_i[get_mod_i[j]], get_t_i[get_mod_i[j]]]
      trad.j[j] = trad.ct[get_c_i[get_mod_i[j]], get_t_i[get_mod_i[j]]]
      modern[get_mod_i[j]] ~ dnorm(modern.j[j], prec)
      trad[get_mod_i[j]] ~ dnorm(trad.j[j], prec)
    }
   for (i in 1:n_unmet){
        unmet[get_unmet_i[i]] ~ dnorm(unmet.ct[get_c_i[get_unmet_i[i]], get_t_i[get_unmet_i[i]]], prec)
   }
   for (k in 1:n_ptot){
    ptot.k[k] = mod.ct[get_c_i[get_ptot_i[k]], get_t_i[get_ptot_i[k]]] + trad.ct[get_c_i[get_ptot_i[k]], get_t_i[get_ptot_i[k]]]
    logit.ptot[k] ~ dnorm(logit(ptot.k[k]), prec)
  }
   
  ",sep="",append=TRUE, file = "model.txt", fill = TRUE)
    
  } else if (!nulldata) {
    # need to check if any updates in bounds min/max used here in Mark's most recent version
    # to decide if still to add periods?
    #  for (j in 1:n_mod){
    #  for (h in 1:getperiod.j[j]) {
    #    trad.jh[j, h] = p.ci[getc.j[j], getest.jf[j, h]] * (1 - R.ci[getc.j[j], getis.jf[j, h]])
    #      modern.jh[j, h] = p.ci[getc.j[j], getest.jf[j, h]] * R.ci[getc.j[j], getis.jf[j, h]]
    #      unmet.jh[j, h] = (1 - p.ci[getc.j[j], getest.jf[j, h]])*(1 / (1 + neg.explogitZ.ci[getc.j[j], getis.jf[j, h]]))
    #}
    #trad.j[j] = 1 / period.j[j] * inprod(trad.jh[j, 1:getperiod.j[j]], partialtime.xj[1:getperiod.j[j], j])
    #modern.j[j] = 1 / period.j[j] * inprod(modern.jh[j, 1:getperiod.j[j]], partialtime.xj[1:getperiod.j[j], j])
    #unmet.j[j] = 1 / period.j[j] * inprod(unmet.jh[j, 1:getperiod.j[j]], partialtime.xj[1:getperiod.j[j], j])
    
    cat("
####
# dms
  for (j in 1:n_mod){
# get_mod_i refers to indices with modern+trad use
   ratios.trad.modern.jn[get_mod_i[j],1:2] ~ dmnorm(mu.jn[get_mod_i[j], ],InvSigma[j,,]) #T.j[j,,])
   InvSigma[j,1:2,1:2] <- inverse(Sigma[j,1:2,1:2])
   Sigma[j,1,2] <- cor.trad.modern.s[get_mod_i[j]]*sqrt(Sigma[j,1,1]*Sigma[j,2,2])
   Sigma[j,2,1] <- cor.trad.modern.s[get_mod_i[j]]*sqrt(Sigma[j,1,1]*Sigma[j,2,2])
   Sigma[j,1,1] <- pow(se_log_r_traditional_no_use[get_mod_i[j]], 2) + pow(nonsample.se.trad.s[get_mod_i[j]],2)
   Sigma[j,2,2] <- pow(se_log_r_modern_no_use[get_mod_i[j]], 2) + pow(nonsample.se.modern.s[get_mod_i[j]],2)
  }
for (i in 1:n_unmet){
# get_unmet_i refers to indices with unmet
 logitratio.yunmet.j[get_unmet_i[i]] ~ dnorm(
   logitratio.yunmet.hat.j[get_unmet_i[i]], 1/(pow(nonsample.se.unmet.s[get_unmet_i[i]],2)+ pow(se_log_r_unmet_no_need[get_unmet_i[i]],2)) )
}
for (k in 1:n_ptot){
    logit.ptot[get_ptot_i[k]] ~ dnorm(logit.ptothat.j[get_ptot_i[k]], tau.sourcetot)
}
tau.sourcetot <- pow(sigma.sourcetot,-2)
# end dms for unmet and trad+modern
# this part refer to union of indices in unmet and modern

  for (j in 1:N){
 modern.j[j] <- mod.ct[get_c_i[j], get_t_i[j]]
 trad.j[j] <- trad.ct[get_c_i[j], get_t_i[j]]
 unmet.j[j] <- unmet.ct[get_c_i[j], get_t_i[j]]
 mu.jn[j,1] <- log(max(0.0000001, q.ij[1,j])/none.adj.j[j])
 mu.jn[j,2] <- log(max(0.0000001, q.ij[2,j])/none.adj.j[j])
 logitratio.yunmet.hat.j[j] <- logit(max(0.0000001,q.ij[3,j])/none.adj.j[j])
 logit.ptothat.j[j] <- logit(max(0.0000001, 1-none.adj.j[j]))
    sump.j[j] <- (trad.j[j]*Vtrad.j[j] + modern.j[j]* Vmodern.j[j]
                        + (1- trad.j[j] - modern.j[j]))
      # old order, 1 is trad!!!
    p.perturb.ij[1,j] <-  trad.j[j]*Vtrad.j[j]/sump.j[j]
      p.perturb.ij[2,j] <-  modern.j[j]* Vmodern.j[j]/sump.j[j]
      p.perturb.ij[3,j] <- unmet.j[j]/sump.j[j]
      p.perturb.ij[4,j] <- (1- trad.j[j] - modern.j[j] - unmet.j[j])/sump.j[j]
      ###Biases
      ##Inclusion of folk methods
      folkbias.j[j] <- step(folk.ind[j]-0.5)*v.folk* p.perturb.ij[3,j]
      ##Absence of probing
      micsbias.j[j] <- step(source.MICS.ind[j]-0.5)* v.mics * p.perturb.ij[1,j]
      ##Sterilization
      modposbias.j[j] <- step(mpos.ind[j]-0.5)*v.mpos* p.perturb.ij[4,j]
      modnegbias.j[j] <- step(mneg.ind[j]-0.5)*v.mneg * p.perturb.ij[2,j]
      ####Perturbed proportions adjusted for biases (1-4)
      q.ij[1,j] <- p.perturb.ij[1,j] - micsbias.j[j] + folkbias.j[j]
      q.ij[2,j] <- p.perturb.ij[2,j] + modposbias.j[j] - modnegbias.j[j]
      q.ij[3,j] <- p.perturb.ij[3,j] + micsbias.j[j] - folkbias.j[j]
      q.ij[4,j] <- p.perturb.ij[4,j] - modposbias.j[j] + modnegbias.j[j]
      none.adj.j[j] <- max(0.0000001, q.ij[3,j] + q.ij[4,j]) #la 2019/3/13
      Vtrad.j[j] <-  (
      V.geo.12i[1,geo.ind[j]]
      * V.age.12i[1,age.ind[j]]
      * V.hw.12i[1,hw.ind[j]]
      * V.emal.12i[1,emal.ind[j]]
      * V.sa.12i[1,sa.ind[j]]
      * V.posbias.12i[1,posbias.ind[j]]
      * V.posage.12i[1, posage.ind[j]]
      * V.negage.12i[1, negage.ind[j]]
      )
      Vmodern.j[j] <- (
      V.geo.12i[2,geo.ind[j]]  ##geographical region
      * V.age.12i[2,age.ind[j]] #Age group different from base (bias unknown)
      * V.hw.12i[2,hw.ind[j]] ##Husband and wives or both
      * V.emal.12i[2,emal.ind[j]] ##Ever married, all women
      * V.sa.12i[2,sa.ind[j]] ##All sexually active
      * V.posbias.12i[2,posbias.ind[j]] ## Non-pregnant/fertile/married SA women
      * V.posage.12i[2, posage.ind[j]] ##Age group with positive bias
      * V.negage.12i[2, negage.ind[j]] ##Age group with negative bias
      )
  }
# add dummy column in case ncol(V) = 1 => V becomes vector!
V.geo.12i[1,max(geo.ind)+1] <- 0
      V.age.12i[1,max(age.ind)+1] <- 0
      V.hw.12i[1,max(hw.ind)+1] <- 0
      V.emal.12i[1,max(emal.ind)+1] <- 0
      V.sa.12i[1,max(sa.ind)+1] <- 0
      V.posbias.12i[1,max(posbias.ind)+1] <- 0
      V.posage.12i[1,max(posage.ind)+1] <- 0
      V.negage.12i[1,max(negage.ind)+1] <- 0
      V.geo.12i[2,max(geo.ind)+1] <- 0
      V.age.12i[2,max(age.ind)+1] <- 0
      V.hw.12i[2,max(hw.ind)+1] <- 0
      V.emal.12i[2,max(emal.ind)+1] <- 0
      V.sa.12i[2,max(sa.ind)+1] <- 0
      V.posbias.12i[2,max(posbias.ind)+1] <- 0
      V.posage.12i[2,max(posage.ind)+1] <- 0
      V.negage.12i[2,max(negage.ind)+1] <- 0
# Multipliers V in [0,inf): geo, emal, hw, age other, sa (for trad only)
####All of these should get a log normal distribution....
        V.sa.12i[1,1] <- 1  ##Sexually active women (trad)
        V.geo.12i[1,1] <- 1 ##Geographical region (trad)
        V.geo.12i[2,1] <- 1 ##Geographical region (mod)
        V.age.12i[1,1] <- 1 ##Age different (trad)
        V.age.12i[2,1] <- 1 ##Age different (mod)
        V.hw.12i[1,1] <- 1  ##Husband/wives (trad)
        V.hw.12i[2,1] <- 1  ##Husband/wives (mod)
        V.emal.12i[1,1] <- 1 ##ever married/ all women (trad)
        V.emal.12i[2,1] <- 1 ##evr married/ all women (mod)
        
for (m in 1:2){
        for (i in 2:ncat.geo){
        V.geo.12i[m,i] ~ dlnorm(0, tau.geo.m[m])
        }
        for (i in 2:ncat.age){
        V.age.12i[m,i] ~ dlnorm(0, tau.geo.m[m])
        }
        tau.geo.m[m] <- pow(sigma.geo.m[m], -2)
        }
        for (i in 2:ncat.sa){
        V.sa.12i[1,i] ~ dlnorm(0, tau.geo.m[1])
        }
        V.sa.12i[2,1] <- 1
        V.posbias.12i[1,1] <- 1
        V.posbias.12i[2,1] <- 1
        V.posage.12i[1,1] <- 1
        V.posage.12i[2,1] <- 1
        V.negage.12i[1,1] <- 1
        V.negage.12i[2,1] <- 1
        # m = 2:
        for (i in 2:ncat.sa){
        W.sa.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
        V.sa.12i[2,i] <- 1+W.sa.12i[2,i]
        }
        for (i in 2:ncat.posbias){
        W.posbias.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
        V.posbias.12i[2,i] <- 1+W.posbias.12i[2,i]
        }
        for (i in 2:ncat.posage){
        W.posage.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
        V.posage.12i[2,i] <-1+W.posage.12i[2,i]
        }
        for (i in 2:ncat.negage){
        W.negage.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
        V.negage.12i[2,i] <- 1/(1+W.negage.12i[2,i])
        }
        tau.pos <- pow(sigma.pos, -2)
        # m=1
        # note: could simplify code and throw out these V's
        for (i in 2:ncat.posbias){
        V.posbias.12i[1,i] <- 1+exp(mu.pos.m[1])
        }
        for (i in 2:ncat.posage){
        V.posage.12i[1,i] <- 1+exp(mu.pos.m[1])
        }
        for (i in 2:ncat.negage){
        V.negage.12i[1,i] <- 1/(1+exp(mu.pos.m[1]))
        }
  ",sep="",append=TRUE, file = "model.txt", fill = TRUE)
  }
if (is_in_union == "Y" & !nulldata) {
  cat("
        for(k in 1:n.training.modonly) {
            logit.ymodonly.j[getj.training.modonly.k[k]] ~ dnorm(logit(q.ij[2,getj.training.modonly.k[k]]), tau.sourcemodonly)
        }
        ",sep="",append=TRUE, file = "model.txt", fill = TRUE)
}
if (is_in_union == "Y") {
  cat("
        for (m in 1:2){
        for (i in 2:ncat.emal){
        V.emal.12i[m,i] ~ dlnorm(0, tau.geo.m[m])
        }
        for (i in 2:ncat.hw){
        V.hw.12i[m,i] ~ dlnorm(0, tau.geo.m[m])
        }
        }
  
  
  ",sep="",append=TRUE, file = "model.txt", fill = TRUE)
} else {
  cat("
  ## m = 1
  for (i in 2:ncat.emal){
    V.emal.12i[1,i] <- 1/(1+exp(mu.pos.m[1]))
  }
  
  ## m = 2
  for (i in 2:ncat.emal){
    W.emal.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
    V.emal.12i[2,i] <- 1/(1+W.emal.12i[2,i])
  }
  
  ## m = 1
  for (i in 2:ncat.hw){
    V.hw.12i[1,i] <- 1+exp(mu.pos.m[1])
  }
  
  ## m = 2
  for (i in 2:ncat.hw){
    W.hw.12i[2,i] ~ dlnorm(mu.pos.m[2], tau.pos)
    V.hw.12i[2,i] <- 1+W.hw.12i[2,i]
  }
  ",sep="",append=TRUE, file = "model.txt", fill = TRUE)
  
}


if (include_ss_data) {
  cat("
      for (k in 1:(n_ss-1)) {
        delta_emu[k] ~ dnorm(
        delta_modern[k], tau_emu[k])
        delta_modern[k] <- mod.ct[1,get_t_ss[k+1]] - mod.ct[1,get_t_ss[k]]
        tau_emu[k] <- pow(se_emu[k], -2)
}
",sep="",append=TRUE, file = "model.txt", fill = TRUE)
}
cat("} # end model",sep="",append=TRUE, file = "model.txt", fill = TRUE)
} # end write model function