Search code examples
rloopslapplystatistics-bootstrap

Looping through function arguments (series of contrasts with multcomp::glht)


I wish to write a function that runs contrasts over a regression model and bootstraps those results to get confidence intervals, looping that function over a list of contrasts.

I have tried for loops nested within functions, lapply, map ... none seem to get me what I want (returns results for either only the first contrast in the list or the last).

For a single contrast from the list of contrasts, the code looks like this:

df <- data.frame(

  H0013301_new_data = c(0,2,3,6,0,4,2,4,8,1),
  drink_stat94_KEYES_2 = c("Heavy","Abstainer","Occasional","Moderate","Abstainer","Occasional","Heavy","Moderate","Moderate","Abstainer"),
  drink_stat02_KEYES_2 = c("Heavy","Abstainer","Occasional","Abstainer","Abstainer","Heavy","Heavy","Moderate","Moderate","Abstainer"),
  drink_stat06_KEYES_2 = c("Occasional","Abstainer","Occasional","Abstainer","Occasional","Heavy","Heavy","Moderate","Moderate","Heavy"),
FIN_weight_survPS_trimmed=
c(.5,2.4,.6,4.8,1.2,.08,.34,.56,1.6,.27)
)

#reordering factors
df$drink_stat94_KEYES_2<-fct_relevel(df$drink_stat94_KEYES_2, "Abstainer", "Occasional", "Moderate", "Heavy")
contrasts(df$drink_stat94_KEYES_2)<-contr.treatment(4,base=1)
df$drink_stat02_KEYES_2<-fct_relevel(df$drink_stat02_KEYES_2, "Abstainer", "Occasional", "Moderate", "Heavy")
contrasts(df$drink_stat02_KEYES_2)<-contr.treatment(4,base=1)
df$drink_stat06_KEYES_2<-fct_relevel(df$drink_stat06_KEYES_2, "Abstainer", "Occasional", "Moderate", "Heavy")
contrasts(df$drink_stat06_KEYES_2)<-contr.treatment(4,base=1)



#defining contrast
c1 <- rbind("A,A,A"=c(1,0,0,0,0,0,0,0,0,0)
                ) 

#defining function to feed to boostrap
fc_2<-function(d,i){
 TrialOutcomeModel_M<-lm(H0013301_new_data ~ drink_stat94_KEYES_2 + drink_stat02_KEYES_2 + drink_stat06_KEYES_2, weights=FIN_weight_survPS_trimmed, data = d[i,]) 
 test <- multcomp::glht(TrialOutcomeModel_M, linfct=c1) 
 return(coef(test))
}
boot_out<-boot(data=df, fc_2, R=500) 
boot.ci(boot_out, type="perc")

But let's assume that instead of just c1, I want to run my function (and boostrap the results) over the following list of contrasts:

c1 <- rbind("A,A,A"=c(1,0,0,0,0,0,0,0,0,0)
                ) 
c2 <- rbind("A,A,O"=c(1,0,0,0,0,0,0,1,0,0)
                ) 
c3 <- rbind("A,A,M"=c(1,0,0,0,0,0,0,0,1,0)
                ) 

c_vector<-list(c1,c2,c3)

Any suggested code for how I would go about this? (P.S. I know that the linfct argument can take a matrix of contrasts, but I'm specifically looking for a loop/lapply solution).


Solution

  • (In the following I'll reference the objects you create in the example code)

    The plan has 2 steps:

    1. preparing a function fun_boot() that takes a contrast object (like c1), and returns a boot object based on it and the df data;

    2. applying that function to the list c_vector of contrasts.

    Consequently, the implementation has 2 elements:

    # [!] Assume all required libraries loaded
    # [!] Assume all necessary data exists
    
    # Step 1
    fun_boot <- function(contrast)
    {
        # Make statistic function
        fun_statistic <- function(d, i)
        {
            TrialOutcomeModel_M <- lm(
               formula = H0013301_new_data ~ drink_stat94_KEYES_2 + drink_stat02_KEYES_2 + drink_stat06_KEYES_2,
               data    = d[i,],
               weights = FIN_weight_survPS_trimmed
            ) 
            test <- multcomp::glht(
               TrialOutcomeModel_M,
               linfct = contrast
            )
            return(coef(test))
        }
    
        # Make boot call (hehe)
        return (boot(
           data      = df,
           statistic = fun_statistic,
           R         = 500
        ))
    }
    
    # Step 2
    boot_out_vector <- lapply(
       X   = c_vector,
       FUN = fun_boot
    )