Implementation of Equations in S+




The following S+ functions use the S+ function ms() to solve for the parameter values explained in the equations section.

bwt and infd are the birth weight and infant mortality variables. It is assumed that there are no missing values or other problematic values and that infd is coded as 0 = alive and 1 = dead.

start.value is a variable which contains the starting values needed for the ms() routine. In the case of the univariate model with no additional covariates, there are eleven parameters.

mix.trunc is a value which indicates where the birth weight data is should be truncated. We often use 500 g.

The Uni.Bwt.Find.PE.Func returns a list of the parameters, the maximum likelihood value, and the convergence flag from ms().

Uni.Bwt.Find.PE.Func <- function(bwt, infd, start.value, mix.trunc`)
{
            Cohort <- cbind(bwt, infd)
            if (mix.trunc.value>0)
                        {
                        bwt <- cohort[, 1]
                        cohort <- cohort[bwt>=mix.trunc.value, ]
                        }
            cat(dim(cohort), "\n\n")
            n1 <- dim(cohort)[1]
            bwt <- cohort[, 1]
            infd <- cohort[, 2]

            joint.PE <- rep(1, 11)
            res.joint <- ms(~ joint.mle.func(bwt, infd, mix.trunc, joint.PE), trace = F, start = list(joint.PE = start.value), control = list(max.step = 1, maxiter = 500, maxfcalls = 2000))
            cat(res.joint$parameters, "\n\n")           

            result <- c(n0, n1, res.joint$parameters, res.joint$value, min.cvg.msg(res.joint$flags[1]))
            names(result) <- c("n0", "n1”, "p.sec", "m.sec", "sd.sec", "m.pri", "sd.pri","a.sec", "b.sec", "c.sec", "a.pri", "b.pri", "c.pri", "joint.mlv", "joint.conv.flag")
            result
}

 

joint.mle.func <-
function(bwt, infd, mix.trunc, joint.PE)
{
            subpop.prob <- mix.cal.subpop.ratio.func(bwt, mix.trunc, joint.PE[1:5])
            grp.ratio <- subpop.prob[, 3]
            bsq <- bwt*bwt
            numerator <- exp(infd * (joint.PE[6] + joint.PE[7] * bwt + joint.PE[8] * bsq))
            denominator <- 1 + exp(joint.PE[6] + joint.PE[7] * bwt + joint.PE[8] * bsq)
            sec.prob <- numerator/denominator * grp.ratio
            sec.prob[numerator == Inf & denominator == Inf] <- grp.ratio[numerator == Inf & denominator == Inf]
            numerator <- exp(infd * (joint.PE[9] + joint.PE[10] * bwt + joint.PE[11] * bsq))
            denominator <- 1 + exp(joint.PE[9] + joint.PE[10] * bwt + joint.PE[11] * bsq)
            pri.prob <- numerator/denominator * (1-grp.ratio)
            pri.prob[numerator == Inf & denominator == Inf] <- 1 - grp.ratio[numerator == Inf & denominator == Inf]
            tot.prob <-  - log(sec.prob + pri.prob) + ( - log(subpop.prob[, 1] + subpop.prob[, 2]))           
            if(length(na.omit(tot.prob)) < length(tot.prob)) cat(joint.PE)
            tot.prob

}