R/fimle.lnl.R

In iMediate: Likelihood Methods for Mediation Analysis

#' Full Information Maximum Likelihood Estimates in Linear M-model and Linear
#' Y-model
#' 
#' Parameter estimates in system of correlated linear M-model and linear
#' Y-model with treatment-mediator interaction using the full information
#' maximum likelihood method.
#' 
#' %% ~~ If necessary, more details than the description above ~~ 
#' P-values are computed from normal distribution.
#' 
#' @param fit.M a fitted model object for mediator. It must be an object
#' generated by function ``\code{lm}''
#' @param fit.Y a fitted model object for outcome. It must be an object
#' generated by function ``\code{lm}''. It can contain treatment-mediator
#' interaction
#' @param treatment a character string of the name of the treatment variable. This variable takes numerical values
#' @param rho a numerical variable specifying the correlation coefficient
#' between the residual of the M-model and the residual of the Y-model. Its
#' range is between \code{-1} and \code{1}
#' @param t0 a reference value for the treatment
#' @param t1 another value for the treatment 
#' @param m a value specifying the level of the mediator. Used for CDE computation 
#' between the residual of the M-model and the residual of the Y-model. Its
#' range is between \code{-1} and \code{1}
#' @return A list containing the following components: \item{M.model}{a data
#' frame containing the results for the M-model } \item{Y.model}{a data frame
#' containing the results for the Y-model } \item{Effects}{a data frame containing 
#' estimated ACME, ADE, Total Effect, and CDE for treatment values \code{t1} and \code{t0}}
#' \item{Variance}{a matrix of 
#' variances and covariances of the parameters estimates }
#' @author Kai Wang \code{<kai-wang@@uiowa.edu>}
#' @references Wang, K. (2019) Maximum likelihood analysis of mediation models with treatment-mediator interaction.
#' Revision submitted.
#' @keywords estimates
#' @examples
#' data("jobs", package = "mediation")
#' 
#' fit.M <- lm(job_seek ~ treat + econ_hard + sex + age, data=jobs)
#' fit.Y <- lm(depress2 ~ treat + job_seek + econ_hard + sex + age, data=jobs)
#' fimle.lnl(fit.M, fit.Y, "treat", rho=0.2)
#' 
#' fit.M <- lm(job_seek ~ treat + econ_hard + sex + age , data=jobs)
#' fit.Y <- lm(depress2 ~ treat*job_seek+ econ_hard + sex + age , data=jobs)
#' fimle.lnl(fit.M, fit.Y, "treat", rho=0.5)
#'
#' @export fimle.lnl
#' @import stats
fimle.lnl = function(fit.M, fit.Y, treatment, rho=0, t0=0, t1=1, m=1) {
     M.name = all.vars(formula(fit.M))[1]
     Y.factors = attr(terms(formula(fit.Y)), "term.labels")
     n = length(residuals(fit.M))
 
     s1.2 = mean(residuals(fit.M)^2)
     s2.2 = mean(residuals(fit.Y)^2)/(1-rho^2)              
     tau = rho*sqrt(s2.2/s1.2)

     M.c = coef(summary(fit.M))
     Y.c = coef(summary(fit.Y))
     for (var in row.names(M.c)){
        	Y.c[var,"Estimate"] = Y.c[var,"Estimate"] + M.c[var,"Estimate"]*tau
     }
     Y.c[M.name,"Estimate"] = Y.c[M.name,"Estimate"] - tau
     L1 = model.matrix(fit.M)
     L2 = model.matrix(fit.Y)

     A = rbind(cbind(crossprod(L1)/s1.2, -rho/sqrt(s1.2*s2.2)*crossprod(L1, L2)), 
               cbind(-rho/sqrt(s1.2*s2.2)*crossprod(L2, L1), crossprod(L2)/s2.2))/(1-rho^2)
     B = matrix(0, nrow = ncol(L1), ncol=2)
     B2 = kronecker(matrix(c(-1/s1.2, 1/s2.2), ncol=2)*rho/2/(1-rho^2)/sqrt(s1.2*s2.2), crossprod(L2, residuals(fit.M))) 
     B = rbind(B, B2)
     C = n/4/(1-rho^2)*matrix(c((2-rho^2)/s1.2^2, -rho^2/s1.2/s2.2, -rho^2/s1.2/s2.2, (2-rho^2)/s2.2^2), ncol=2)

     VC = solve(A - B %*% solve(C, t(B)))
     SEs = sqrt(diag(VC))
     M.SEs = SEs[1:nrow(M.c)]
     Y.SEs = SEs[-(1:nrow(M.c))]
     for (var in row.names(M.c)) M.c[var,"Std. Error"] = M.SEs[names(M.SEs) == var]
     for (var in row.names(Y.c)) Y.c[var,"Std. Error"] = Y.SEs[names(Y.SEs) == var]
 
     M.c[, "t value"] = M.c[,"Estimate"]/M.c[,"Std. Error"]
     M.c[, "Pr(>|t|)"] = signif(2*pnorm(abs(M.c[,"t value"]), lower.tail=FALSE), 4)
     Y.c[, "t value"] = Y.c[,"Estimate"]/Y.c[,"Std. Error"]
     Y.c[, "Pr(>|t|)"] = signif(2*pnorm(abs(Y.c[,"t value"]), lower.tail=FALSE), 4)

    colnames(M.c) = colnames(Y.c) = c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
    tt1 = paste(treatment, M.name, sep=":")
    tt2 = paste(M.name, treatment, sep=":")
    kap = 0
    inter.term = FALSE
    if ((tt1 %in% Y.factors) | (tt2 %in% Y.factors)) {
        if (tt1 %in% Y.factors) TM = tt1
    	    else TM = tt2
        kap = Y.c[TM, "Estimate"]
        inter.term = TRUE
    	}
    
    delta = function(tt){
    	    M.c[treatment, "Estimate"]*(Y.c[M.name, "Estimate"] + kap*tt)*(t1-t0)
    }
    
    zeta = function(tt){
        tmp = mean(predict(fit.M)) + M.c[treatment, "Estimate"]*(tt - mean(L1[, treatment]))
        (Y.c[treatment, "Estimate"] + kap*tmp)*(t1-t0)
    }

    total = delta(t1) + zeta(t0)
    cde = (Y.c[treatment, "Estimate"] + kap*m)*(t1-t0)
    
    var.u = function(tt){
    	    	u.M = rep(0, dim(L1)[2])
    	    	names(u.M) = dimnames(L1)[[2]]
    	    	u.M[treatment] = (Y.c[M.name, "Estimate"]+kap*tt)*(t1-t0)
    	    	u.Y = rep(0, dim(L2)[2])
    	    	names(u.Y) = dimnames(L2)[[2]]
    	    	u.Y[M.name] = M.c[treatment, "Estimate"]*(t1-t0)
        if (inter.term) u.Y[TM] = M.c[treatment, "Estimate"]*tt*(t1-t0)
        uu = c(u.M, u.Y)
        crossprod(uu, VC %*% uu)
     }
    	    
    Xbar = colMeans(L1[, !(dimnames(L1)[[2]] %in% c("(Intercept)", treatment))]) # vector of covariate means
    var.v = function(tt){
    	    v.M = kap*c(1, tt, Xbar)*(t1-t0)
    	    v.Y = rep(0, dim(L2)[2])
    	    names(v.Y) = dimnames(L2)[[2]]
    	    v.Y[treatment] = 1*(t1-t0)
    	    tmp = mean(predict(fit.M)) + M.c[treatment, "Estimate"]*(tt - mean(L1[, treatment]))
        if (inter.term) v.Y[TM] = tmp*(t1-t0)
        vv = c(v.M, v.Y)
        crossprod(vv, VC %*% vv)
     }

    w.M = c(kap, Y.c[M.name, "Estimate"] + kap*(t1+t0), Xbar)*(t1-t0)
   	w.Y = rep(0, dim(L2)[2])
   	names(w.Y) = dimnames(L2)[[2]]
   	w.Y[treatment] = 1*(t1-t0)
   	w.Y[M.name] = M.c[treatment, "Estimate"]*(t1-t0)
    if (inter.term) w.Y[TM] = (M.c["(Intercept)", "Estimate"] + M.c[treatment, "Estimate"]*(t1+t0))*(t1-t0)
    ww = c(w.M, w.Y)
    var.w = crossprod(ww, VC %*% ww)

   	s.M = rep(0, dim(L1)[2])
   	s.Y = rep(0, dim(L2)[2])
   	names(s.Y) = dimnames(L2)[[2]]
   	s.Y[M.name] = 1*(t1-t0)
    if (inter.term) s.Y[TM] = m*(t1-t0)
    ss = c(s.M, s.Y)
    var.cde = crossprod(ss, VC %*% ss)

    tmp = c(delta(t0), delta(t1), zeta(t0), zeta(t1), total, cde)
    SE = sqrt(c(var.u(t0), var.u(t1), var.v(t0), var.v(t1), var.w, var.cde))
    effects = data.frame(Effect = tmp, SE = SE, 
                         LB = tmp - qnorm(0.975)*SE, 
                         UB = tmp + qnorm(0.975)*SE)
    effects$zvalue = effects$Effect/effects$SE
    effects$pvalue = signif(2*pnorm(abs(effects$zvalue), lower.tail=FALSE), 4)

    names(effects) = c("Effect", "SE", "95% LB", "95% UB", "z value", "Pr(>|z|)")
    row.names(effects) = c(paste("ACME (t=", t0, ")", sep=""),
                           paste("ACME (t=", t1, ")", sep=""),
                           paste("ADE (t=", t0, ")", sep=""), 
                           paste("ADE (t=", t1, ")", sep=""), 
                           "Total Effect",
                           paste("CDE (m=", m, ")", sep=""))
 
     list(M.model = as.data.frame(M.c), Y.model = as.data.frame(Y.c), Effects = effects, Variance = VC)
}