R/larf.R

In LARF: Local Average Response Functions for Instrumental Variable Estimation of Treatment Effects

#######################################
### LARF: Local Response Functions in R
#######################################

##############################
### The high-interface
##############################

larf <- function(formula, treatment = NULL, instrument = NULL, data, method = "LS", AME = FALSE, optimizer = "Nelder-Mead", zProb = NULL ) {
  
  ## set up a model.frame
  if(missing(data)) data <- environment(formula)
  
  if (!is.null(treatment)) {
    mf <- model.frame(formula = formula, data = data)
    ## extract response, covaraites, treatment, and instrument
    Y <- model.response(mf)
    X <- model.matrix(attr(mf, "terms"), data = mf)
    D <- treatment
    Z <- instrument
  }
  else {  
    f <- Formula::Formula(formula)
    formula <- f
    Y <- Formula::model.part(f, data = data, lhs = 1)[,1]
    D <- Formula::model.part(f, data = data, lhs = 1)[,2]
    X <- model.matrix(f, data = data, rhs = 1)
    Z <- model.matrix(f, data = data, rhs = 2)[,-1]       
  }
  
  ## call default interface
  est <- larf.fit(Y, X, D, Z, method, AME, optimizer, zProb) 

  est$call <-match.call()
  est$formula <- formula

  return(est)
}


##############################
### Function to implement LARF
##############################

larf.fit <- function(Y, X, D, Z, method, AME, optimizer, zProb) {
  y <- as.matrix(Y)
  X <- as.matrix(X)
  d <- as.matrix(D)
  z <- as.matrix(Z)
  
  n <- length(y)   
  
  outcome <- floor((colSums(y==1 | y==0))/n) 

  # Get the weights kappa
  if ( is.null(zProb) ) {
  eqA <- glm(z ~ X - 1, family=binomial(probit))
  gamma <- as.matrix(summary(eqA)$coefficients[,1])
  Phi <- eqA$fitted
  } else {
    Phi <- zProb
    lv <- Z*(1-D)/Phi^2 - D*(1-Z)/(1-Phi)^2
  }
  
  nc_ratio <- rep(0,n)
  nc_ratio[d==0 & z==1] <- 1/Phi[d==0 & z==1]
  nc_ratio[d==1 & z==0] <- 1/(1-Phi[d==1 & z==0])
  kappa <- 1 - nc_ratio
    
  #  LARF for a binary outcome
  if(outcome == 1) {
    
    # Get initial s of the parameters
    eqB <- glm(y ~ d + X - 1, family=binomial(probit))
    theta <- as.matrix(summary(eqB)$coefficients[,1])
    b1 <- theta
    
    ## Least squares (Abadie 2003: equation 8)
    if (method == "LS") {
      step <- 1
      gtol <- 0.0001
      dg <- 1
      k <- 1
      while( dg > gtol & k <= 1000){
        b0 <- b1
        u <- y - pnorm(cbind(d,X)%*%b0) # The residual
        g <- t(cbind(d,X)*as.vector(kappa*dnorm(cbind(d,X)%*%b0)))
        delta <- step*solve(tcrossprod(g))%*%g%*%u # The Gauss-Newton Method
        b1 <- b0 + delta
        dg <- crossprod(delta,g)%*%crossprod(g,delta)
        k <- k + 1
      }
      b <- b1    
  
      if ( is.null(zProb) ) {   # Get the SE using theorem 4.2
        u <- y - pnorm(cbind(d,X)%*%b)
        
        dM_theta <- (cbind(d,X)%*%b)*dnorm(cbind(d,X)%*%b)*u+(dnorm(cbind(d,X)%*%b))^2
        M_theta <- t(cbind(d,X)*as.vector(kappa*dM_theta))%*%cbind(d,X)
        
        derkappa <- rep(0,n)
        derkappa[z==1&d==0] <- 1/(Phi[z==1&d==0])^2
        derkappa[z==0&d==1] <- (-1)/(1-Phi[z==0&d==1])^2
        
        M_gamma <- -t(cbind(d,X)*as.vector(dnorm(cbind(d,X)%*%b)*u*matrix(derkappa)*dnorm(X%*%gamma)))%*%X
        
        lambda <- rep(0,n)
        lambda[z==0] <- (-1)*dnorm(X[z==0,]%*%gamma)/(1-pnorm(X[z==0,]%*%gamma))
        lambda[z==1] <- dnorm(X[z==1,]%*%gamma)/pnorm(X[z==1,]%*%gamma)
        
        H_gamma <- t(X*lambda^2)%*%X
        
        # Get the variance
        S_1 <- t(cbind(d,X)*as.vector((dnorm(cbind(d,X)%*%b)*u*kappa)^2))%*%cbind(d,X)
        S_2 <- M_gamma%*%tcrossprod(solve(H_gamma),(M_gamma))
        S_3 <- t(t(M_gamma%*%tcrossprod(solve(H_gamma),X))*as.vector(lambda*kappa*(-u)*dnorm(cbind(d,X)%*%b)))%*%cbind(d,X)
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21) %*% S %*% solve(M_theta,tol=1e-21)
        se <- as.matrix(sqrt(diag(V)))
      } else {                  # Get the SE using theorem 4.5
        u <- y - pnorm(cbind(d,X)%*%b)      
        dM_theta <- (cbind(d,X)%*%b)*dnorm(cbind(d,X)%*%b)*u+(dnorm(cbind(d,X)%*%b))^2
        M_theta <- t(cbind(d,X)*as.vector(kappa*dM_theta))%*%cbind(d,X)
        
        # Get the variance
        S_1 <- t(cbind(d,X)* as.vector((dnorm(cbind(d,X)%*%b)*u*kappa)^2))%*%cbind(d,X)
        S_2 <- t(cbind(d,X)*as.vector((dnorm(cbind(d,X)%*%b)*u*lv*(Z-Phi))^2))%*%cbind(d,X)
        S_3 <- t(cbind(d,X)*as.vector((dnorm(cbind(d,X)%*%b)*u)^2*kappa*lv*(Z-Phi)))%*%cbind(d,X)
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21) %*% S %*% solve(M_theta,tol=1e-21)
        se <- as.matrix(sqrt(diag(V)))      
      }
    }
    
    ## Maximum Likelihood (Abadie 2003: equation 10)
    if (method == "ML") {
      target<-function(b){
        -sum(kappa*(y*log(pnorm(cbind(d,X)%*%b))+(1-y)*log(pnorm(-cbind(d,X)%*%b))))
      }
      b2<-optim(theta, target, method = optimizer, hessian = FALSE, control=list(maxit=1e8,abstol=0.1e-8))$par
      b <- b2
      
      if ( is.null(zProb) ) {   # Get the SE using theorem 4.2
        u <- pnorm(cbind(d,X)%*%b)
        r <- dnorm(cbind(d,X)%*%b)
        dM_theta <- -r*(y*(r+cbind(d,X)%*%b*u)/u^2 + (1-y)*(r-cbind(d,X)%*%b*(1-u))/(1-u^2))
        M_theta <- t(cbind(d,X)*as.vector(kappa*dM_theta))%*%cbind(d,X)
        
        derkappa <- rep(0,n)
        derkappa[z==1&d==0] <- 1/(Phi[z==1&d==0])^2
        derkappa[z==0&d==1] <- (-1)/(1-Phi[z==0&d==1])^2
        
        M_gamma <- t(cbind(d,X)*as.vector((y*r/u-(1-y)*r/(1-u))*matrix(derkappa)*dnorm(X%*%gamma)))%*%X
        
        lambda <- rep(0,n)
        lambda[z==0] <- (-1)*dnorm(X[z==0,]%*%gamma)/(1-pnorm(X[z==0,]%*%gamma))
        lambda[z==1] <- dnorm(X[z==1,]%*%gamma)/pnorm(X[z==1,]%*%gamma)
        
        H_gamma <- t(X*lambda^2)%*%X        
      
        S_1 <- t(cbind(d,X)*as.vector(((y*r/u-(1-y)*r/(1-u))*kappa)^2))%*%cbind(d,X)
        S_2 <- M_gamma%*%tcrossprod(solve(H_gamma),(M_gamma))
        S_3 <- t(t(M_gamma%*%tcrossprod(solve(H_gamma),X))*as.vector(lambda*kappa*(y*r/u-(1-y)*r/(1-u))))%*%cbind(d,X)
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21) %*% S %*% solve(M_theta,tol=1e-21)
        se <- as.matrix(sqrt(diag(V)))
      } else {                  # Get the SE using theorem 4.5
        u <- pnorm(cbind(d,X)%*%b)
        r <- dnorm(cbind(d,X)%*%b)
        dM_theta <- -r*(y*(r+cbind(d,X)%*%b*u)/u^2 + (1-y)*(r-cbind(d,X)%*%b*(1-u))/(1-u^2))
        M_theta <- t(cbind(d,X)*as.vector(kappa*dM_theta))%*%cbind(d,X)
        
        S_1 <- t(cbind(d,X)*as.vector(((y*r/u-(1-y)*r/(1-u))*kappa)^2))%*%cbind(d,X)
        S_2 <- t(cbind(d,X)*as.vector(((y*r/u-(1-y)*r/(1-u))*lv*(Z-Phi))^2))%*%cbind(d,X)
        S_3 <- t(cbind(d,X)*as.vector((y*r/u-(1-y)*r/(1-u))^2*kappa*lv*(Z-Phi)))%*%cbind(d,X)
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21) %*% S %*% solve(M_theta,tol=1e-21)
        se <- as.matrix(sqrt(diag(V)))      
      }
    }

    # Marginal effects at means
    cat <- floor((colSums(cbind(d,X)==1 | cbind(d,X)==0))/n) 

    if (AME == 0) {    
      wbar <- apply(kappa*cbind(d,X),2,mean) / mean(kappa)
      db <- as.vector(dnorm(wbar%*%b))*b 
      # P781 in Green
      G <- as.vector(dnorm(wbar%*%b))*(diag(dim(cbind(d,X))[2])-as.vector(wbar%*%b)*b%*%wbar)      
      dV <- G%*%V%*%t(G)
      dse <- as.matrix(sqrt(diag(dV)))  
      
      for(i in 1:length(cat)){
        if (cat[i]==1) { 
          wbar1 <- wbar
          wbar1[i] <- 1
          wbar0 <- wbar
          wbar0[i] <- 0
          db[i] <- pnorm(crossprod(matrix(wbar1),b)) - pnorm(crossprod(matrix(wbar0),b))
          g <- as.vector(dnorm(crossprod(matrix(wbar1),b)))*wbar1 - as.vector(dnorm(crossprod(matrix(wbar0),b)))*wbar0
          dse[i] <- sqrt(g%*%V%*%matrix(g))
          dV[i,i] <- dse[i]^2
        }
      }      
    } 
    
    # Average marginal effects
    if (AME == 1) {
      wbar <- cbind(d,X)     
      db <- mean(dnorm(as.matrix(wbar)%*%b))*b
      P <- dim(cbind(d,X))[2]
      G <- matrix(0, P, P)
      for (i in 1:n) {
        G <- G + as.vector(dnorm(wbar[i,]%*%b))*(diag(P)- as.vector(wbar[i,]%*%b)*b%*%wbar[i,])   
      }
      G <- G / n      
      dV <- G%*%V%*%t(G)
      dse <- as.matrix(sqrt(diag(dV)))
      
      for(i in 1:length(cat)){
        if (cat[i]==1) { 
          wbar1 <- wbar
          wbar1[,i] <- 1
          wbar0 <- wbar
          wbar0[,i] <- 0
          db[i] <- mean(pnorm(as.matrix(wbar1)%*%b) - pnorm(as.matrix(wbar0)%*%b))
          g <- matrix(0, 1, P)
          for (j in 1:n) {
            g <- g + as.vector(dnorm(wbar1[j,]%*%b))*wbar1[j,] - as.vector(dnorm(wbar0[j,]%*%b))*wbar0[j,]
          }
          g <- g / n 
          dse[i] <- sqrt(g%*%V%*%matrix(g))
          dV[i,i] <- dse[i]^2
        }
      }    
    } 
    
    # Report output
    b<-as.matrix(b)
    se<-as.matrix(se)
    tratio<- as.matrix(b/se)
    ptratio <- as.matrix(2 * pt(abs(tratio), df = n-nrow(b), lower.tail=FALSE ))
    db<-as.matrix(db)
    dse<-as.matrix(dse)
    dtratio<- as.matrix(db/dse)
    dptratio <- as.matrix(2 * pt(abs(dtratio), df = n-nrow(b), lower.tail=FALSE ))
    V<-as.matrix(V)
    fitted.values <- pnorm(as.matrix(cbind(D,X))%*%b)
    residuals <- Y - fitted.values
    
    call<-match.call()
    out <- list(call = call, coefficients = b, SE = se, tratio = tratio, ptratio = ptratio, MargEff = db, MargSE = dse, 
           dtratio = dtratio, dptratio = dptratio, vcov = V, fitted.values  = fitted.values, residuals = residuals, outcome = outcome)
    
    rownames(out$coefficients) <- c("Treatment", colnames(X))
    rownames(out$coefficients) <- c(rownames(out$coefficients)[1], gsub("^X", "", rownames(out$coefficients))[-1])
    rownames(out$SE)<-rownames(out$coefficients)
    rownames(out$tratio)<-rownames(out$coefficients)
    rownames(out$ptratio)<-rownames(out$coefficients)
    rownames(out$MargEff)<-rownames(out$coefficients)
    rownames(out$MargSE)<-rownames(out$coefficients)
    rownames(out$vcov)<-rownames(out$coefficients)
    
    class(out)<-"larf"
    return(out)   
  } # End of the binary LARF
  
    
  # For a continuous outcome
  if(outcome != 1) {
    
    if (method == "LS") {     # WLS  of the coefficients
      # b1 <- solve(crossprod(cbind(d,X), DK) %*% cbind(d,X)) %*% crossprod(cbind(d,X), DK) %*% y
      b1 <- solve ( t(cbind(d,X) * kappa) %*% cbind(d,X)) %*% t(cbind(d,X) * kappa)  %*% y
      b <- b1
      
      if ( is.null(zProb) ) {   # Get the SE using theorem 4.2
        u <- y-(cbind(d,X)%*%b)        
        # M_theta <- crossprod(cbind(d,X),diag(as.vector(kappa),n))%*%cbind(d,X)
        M_theta <- t(cbind(d,X) * kappa)%*%cbind(d,X)
        
        derkappa <- replicate(n,0)
        derkappa[z==1&d==0] <- 1/(Phi[z==1&d==0])^2
        derkappa[z==0&d==1] <- (-1)/(1-Phi[z==0&d==1])^2
        
        #M_gamma <- -crossprod(cbind(d,X),diag(as.vector(u*matrix(derkappa)*dnorm(X%*%gamma)),n))%*%X
        M_gamma <- -t( cbind(d,X) * as.vector((u*matrix(derkappa)*dnorm(X%*%gamma))) ) %*% X        
        
        lambda <- replicate(n,0)
        lambda[z==0] <- (-1)*dnorm(X[z==0,]%*%gamma)/(1-pnorm(X[z==0,]%*%gamma))
        lambda[z==1] <- dnorm(X[z==1,]%*%gamma)/pnorm(X[z==1,]%*%gamma)
        #H_gamma <- crossprod(X,diag(lambda^2,n))%*%X
        H_gamma <- t(X *lambda^2) %*% X
        
        # 1st part of "a^2+b^2+2a*b" variance
        #S_1 <- crossprod(cbind(d,X),diag(as.vector((u*kappa)^2),n))%*%cbind(d,X)
        S_1 <- t(cbind(d,X)*as.vector((u*kappa)^2))%*%cbind(d,X)
        
        # 2nd part
        S_2 <- M_gamma%*%tcrossprod(solve(H_gamma),M_gamma)
        
        # 3rd part, AKA, the a*b part
        #S_3 <- M_gamma%*%tcrossprod(solve(H_gamma),X)%*%diag(as.vector(lambda*kappa*(-u)),n)%*%cbind(d,X)
        S_3 <- t(t(M_gamma%*%tcrossprod(solve(H_gamma),X))*as.vector(lambda*kappa*(-u))) %*%cbind(d,X)
        
        # Add them togather
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21)%*%S%*%solve(M_theta,tol=1e-21)
        se <- sqrt(diag(V))
      } else {  # Get the SE using theorem 4.5
        u <- y-(cbind(d,X)%*%b)      
        #M_theta <- crossprod(cbind(d,X),diag(as.vector(kappa),n))%*%cbind(d,X) 
        M_theta <- t(cbind(d,X) * kappa)%*%cbind(d,X) 
        
        #S_1 <- crossprod(cbind(d,X),diag(as.vector((u*kappa)^2),n))%*%cbind(d,X)
        #S_2 <- crossprod(cbind(d,X),diag(as.vector((u*lv*(Z-Phi))^2),n))%*%cbind(d,X)
        #S_3 <- crossprod(cbind(d,X),diag(as.vector(u^2*kappa*lv*(Z-Phi)),n))%*%cbind(d,X)
        
        S_1 <- t(cbind(d,X)*as.vector((u*kappa)^2))%*%cbind(d,X)
        S_2 <- t(cbind(d,X)*as.vector((u*lv*(Z-Phi))^2))%*%cbind(d,X)
        S_3 <- t(cbind(d,X)*as.vector(u^2*kappa*lv*(Z-Phi)))%*%cbind(d,X)
        S <- S_1 + S_2 + S_3 + t(S_3)               
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21)%*%S%*%solve(M_theta,tol=1e-21)
        se <- sqrt(diag(V))            
      }  
    }
    
    if(method=="ML"){      
      targetB<-function(para){        
         sd<-para[length(para)]
         b<-para[1:(length(para)-1)]
         pb <- dnorm(y-cbind(d,X)%*%b,sd)
         pb[pb==0] <- .Machine$double.xmin
         -sum(kappa*log(pb))        
      }
      #DK <- diag(as.vector(kappa), n)
      #theta <- as.vector(solve(crossprod(cbind(d,X), DK) %*% cbind(d,X)) %*% crossprod(cbind(d,X), DK) %*% y)
      theta <- as.vector(solve ( t(cbind(d,X) * kappa) %*% cbind(d,X)) %*% t(cbind(d,X) * kappa)  %*% y)
      
      MLEparameters<-optim(c(theta,sd(y)), targetB, method = optimizer, hessian = FALSE, control=list(maxit=1e8,abstol=0.1e-8))$par
      b2 <- matrix(MLEparameters[1:(length(MLEparameters)-1)])
      b <- b2
      
      # In fact, the sd's can be ignored in the following calculations.
      sd<-MLEparameters[length(MLEparameters)]
      
      if ( is.null(zProb) ) {   # Get the SE using theorem 4.2
        u <- y-cbind(d,X)%*%b
        #M_theta <- -crossprod(cbind(d,X),diag(as.vector(kappa),n))%*%cbind(d,X)/sd^2
        M_theta <- -t(cbind(d,X) * kappa)%*%cbind(d,X)/sd^2
        
        derkappa <- replicate(n,0)
        derkappa[z==1&d==0] <- 1/(Phi[z==1&d==0])^2
        derkappa[z==0&d==1] <- (-1)/(1-Phi[z==0&d==1])^2
        
        #M_gamma <- -crossprod(cbind(d,X),diag(as.vector(u*matrix(derkappa)*dnorm(X%*%gamma)),n))%*%X
        M_gamma <- t(cbind(d,X)*as.vector(u*matrix(derkappa)*dnorm(X%*%gamma)))%*%X/sd^2
        
        lambda <- replicate(n,0)
        lambda[z==0] <- (-1)*dnorm(X[z==0,]%*%gamma)/(1-pnorm(X[z==0,]%*%gamma))
        lambda[z==1] <- dnorm(X[z==1,]%*%gamma)/pnorm(X[z==1,]%*%gamma)
        #H_gamma <- crossprod(X,diag(lambda^2,n))%*%X
        H_gamma <- t(X*lambda^2)%*%X
        
        S_1 <- t(cbind(d,X)*as.vector((u*kappa)^2))%*%cbind(d,X)/sd^4
        S_2 <- M_gamma%*%tcrossprod(solve(H_gamma),M_gamma)
        S_3 <- t(t(M_gamma%*%tcrossprod(solve(H_gamma),X))*as.vector(lambda*kappa*u))%*%cbind(d,X)/sd^2
        S <- S_1 + S_2 + S_3 + t(S_3)
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21)%*%S%*%solve(M_theta,tol=1e-21)
        se <- sqrt(diag(V))         
      }  else{  # Get the SE using theorem 4.5. In fact, the sd's can be ignored.
        u <- y-cbind(d,X)%*%b
        #M_theta <- -crossprod(cbind(d,X),diag(as.vector(kappa),n))%*%cbind(d,X)/sd^2
        M_theta <- -t(cbind(d,X) * kappa)%*%cbind(d,X)/sd^2
        
        S_1 <- t(cbind(d,X)*as.vector((u*kappa)^2))%*%cbind(d,X)/sd^4
        S_2 <- t(cbind(d,X)*as.vector((u*lv*(Z-Phi))^2))%*%cbind(d,X)/sd^4
        S_3 <- t(cbind(d,X)*as.vector(u^2*kappa*lv*(Z-Phi)))%*%cbind(d,X)/sd^4
        S <- S_1 + S_2 + S_3 + t(S_3)               
        
        # The variance and SE
        V <- solve(M_theta,tol=1e-21)%*%S%*%solve(M_theta,tol=1e-21)
        se <- sqrt(diag(V))            
      }  
  }
        
    # Report output
    b<-as.matrix(b)
    se<-as.matrix(se)
    tratio<- as.matrix(b/se)
    ptratio <- as.matrix(2 * pt(abs(tratio), df = n-nrow(b), lower.tail=FALSE ))
    V<-as.matrix(V)
    fitted.values <- as.matrix(cbind(D,X)%*%b)
    residuals <- Y - fitted.values
    
    call<-match.call()
    out <- list(call = call, coefficients = b, SE = se, tratio = tratio, ptratio = ptratio, vcov = V, 
                fitted.values  = fitted.values, residuals = residuals, outcome = outcome)
    
    rownames(out$coefficients) <- c("Treatment", colnames(X))
    rownames(out$SE)<-rownames(out$coefficients)
    rownames(out$tratio)<-rownames(out$coefficients)
    rownames(out$ptratio)<-rownames(out$coefficients)
    rownames(out$vcov)<-rownames(out$coefficients)
    
    class(out)<-"larf"
    return(out)
  } # End of the linear LARF
}
  

##############################
### Generic methods
##############################

print.larf <- function(x, digits = 4, ...)
{  
  cat("Call:\n")
  print(x$call)
  
  est <- cbind(x$coefficients, x$SE)
  colnames(est) <- c("Coefficients", "SE")
  
  cat("\nEstimates:\n")
  print.default(format(est, digits = digits), quote = FALSE)
}


summary.larf <- function(object, ...)
{  
  if (object$outcome == 1) {
    
    TAP<-cbind(Estimate = coef(object),
               SE = object$SE,
               ptratio = object$ptratio,
               MargEff=object$MargEff,
               MargSE=object$MargSE,
               MargP =object$dptratio )
    
    if (is.null(object$call$AME)) object$call$AME <- FALSE      
    if(object$call$AME == 0) { colnames(TAP)<-c("Estimate","SE", "P", "MEM","MEM-SE", "MEM-P") }
    if(object$call$AME == 1) { colnames(TAP)<-c("Estimate","SE", "P", "AME","AME-SE", "AME-P")}
    res <- list(call=object$call, coefficients=TAP)       
  } else {
    TAP<-cbind(Estimate = coef(object), SE = object$SE, ptratio = object$ptratio )
    colnames(TAP)<-c("Estimate","SE", "P")  
    res <- list(call=object$call, coefficients=TAP)  
  }
  class(res) <- "summary.larf"
  return(res)
}


print.summary.larf <- function(x, digits = 4, ...)
{
  cat("Call:\n")
  print(x$call)
  
  cat("\n")
  print.default(round(x$coefficients, digits = digits), quote = FALSE)
}


vcov.larf <- function(object, ...)
{  
  colnames(object$vcov) <- rownames(object$vcov)
  cat("\nCovariance Matrix for the Parameters:\n")
  object$vcov
}


predict.larf <- function(object, newCov, newTreatment, ...) {
  if(object$outcome== 1) {
    object$predicted.values <- pnorm(as.matrix(cbind(newTreatment,newCov))%*%coef(object))
  } 
  else {
    object$predicted.values <- as.matrix(cbind(newTreatment,newCov)%*%coef(object))  
  }
  colnames(object$predicted.values) <- "predicted.values"
  return(object$predicted.values )
}