# This is the "rube" winBugs model for Bayesian fitting of 4-parameter
# logistic trajectories for repeated measures longitudinal data with
# outcomes running strictly between zero and an upper measurement limit.
#
# Random per-subject effects are fit for all four parameters and
# covariates are allowed for the "rate" and "midpoint" parameters.
#  Outliers are accomodated by using an error distribution that follows
# a T rather than a Gaussian distribution.
#
# M is the number of subjects.
# Y is the outcome.  TEST.MAX is its upper limit.
# The per-subject trajectory parameters are yMinFac, yMaxFac,
# tMid, and Rate.  The first two are on a 0-1 scale and represent
# individual curve min and max value when multiplied by TEST.MAX.
# subj[] holds the subject number for each corresponding outcome measurement.
# MID.VARS and RATE.VARS are the covariates for the midpoint and rate
# parameters.
# Other tokens in all upper case are constants or fixed hyperprior
# values.
# Camel case token other than i and j are model parameters.
#
model {
    for (i in 1:M) {
        yMinFac[i] ~ dbeta(minAlpha0, minBeta0)
        yMaxFac[i] ~ dbeta(maxAlpha0, maxBeta0)

        # time of midpoint
        tmptMid[i] <- meantMid + LC(b,MID.VARS,tMid,i)
        tMid[i] ~ dt(tmptMid[i], prectMid, MTDF=3)I(30,180)

        # Rate 
        tmpRate[i] <- meanRate + LC(b,RATE.VARS,Rate,i)
        Rate[i] ~ dt(tmpRate[i], precRate, RTDF=3)I(0,12)
    }

    for (j in 1:NN ) {
        mu[j] <- (TEST.MAX=100)*yMinFac[subj[j]] + 
              ((TEST.MAX=100)*(yMaxFac[subj[j]] - yMinFac[subj[j]]) ) /
            ( 1 + exp( ( X[j] - tMid[subj[j]] ) / Rate[subj[j]] ) ) 
        Y[j] ~ dt( mu[j], precY, TDF=3)
    }
      
    logvarY ~ dunif(LOWLIM.VARY=-2, HILIM.VARY=6)
    precY <- 1/exp(logvarY)
    sdY <- pow(precY, -0.5)

    minAoverB0 ~ dgamma(MIN.ALPHA.P1=2, MIN.ALPHA.P2=1.5)
    minBeta0 ~ dgamma(MIN.BETA.P1=2, MIN.BETA.P2=0.1)
    minAlpha0 <- minAoverB0*minBeta0
    maxAoverB0 ~ dgamma(MAX.ALPHA.P1=2, MAX.ALPHA.P2=0.5)
    maxBeta0 ~ dgamma(MAX.BETA.P1=2, MAX.BETA.P2=0.1)
    maxAlpha0 <- maxAoverB0*maxBeta0

    meantMid ~ dnorm(MEAN.MID=80, PREC.MID=0.01)###I(30,180)
    FOR(b, MID.VARS, tMid, ? ~ dnorm(MEAN.B.MID=0, PREC.B.MID=0.05))

    meanRate ~ dnorm(MEAN.RATE=1, PREC.RATE=0.01)I(0,8)
    FOR(b, RATE.VARS, Rate, ? ~ dnorm(MEAN.B.RATE=0, PREC.B.RATE=0.01))

    logvartMid ~ dunif(LOWLIM.VARMID=-2, HILIM.VARMID=6)
    prectMid <- 1/exp(logvartMid)
    sdtMid <- pow(prectMid, -0.5)
    logvarRate ~ dunif(LOWLIM.VARRATE=-5, HILIM.VARRATE=2)
    precRate <- 1/exp(logvarRate)
    sdRate <- pow(precRate, -0.5)

    # Min/Max
    minMean <- mean(yMinFac[])
    maxMean <- mean(yMaxFac[])
    minSD <- sd(yMinFac[])
    maxSD <- sd(yMaxFac[])
}