using lapply to create t-test table

Question

i want to get t-tests between two populations (in or out of treatment group (1 or 0 in sample data below, respectively)) across a number of variables, and for different studies, all of which are sitting in the same dataframe. In the sample data below, I want to generate t-tests for all variables (in sample data: Age, Dollars, DiseaseCnt) between the 1/0 Treatment group. I want to run these t-tests, by Program, rather than across the population. I have the logic to generate the t-tests. However, I need assistance with the final step of extracting the appropriate parts from the function & creating something easily digestable.

Ultimately, what I want is: a table of t-stats, p-values, variable that t-test was performed on, and program for which variable was tested.

DT<-data.frame(
               Treated=sample(0:1,1000,replace=T)
              ,Program=c('Program A','Program B','Program C','Program D')
              ,Age=as.integer(rnorm(1000,mean=65,sd=15))
              ,Dollars=as.integer(rpois(1000,lambda=1000))
              ,DiseaseCnt=as.integer(rnorm(1000,mean=5,sd=2)) )

progs<-unique(DT$Program) # Pull program names
vars<-names(DT)[3:5] # pull variables to run t tests

test<-lapply(progs, function(i)
          tt<-lapply(vars, function(j) {t.test( DT[DT$Treated==1 & DT$Program == i,names(DT)==j] 
                                                ,DT[DT$Treated==0 & DT$Program == i,names(DT)==j]
                                                ,alternative = 'two.sided'  ) 
              list(j,tt$statistic,tt$p.value)  }
                 ) ) 
  # nested lapply produces results in list format that can be binded, but complete output w/ both lapply's is erroneous

Answer 1

You should convert it into a data.table first. (In my code I call your original table DF):

DT <- as.data.table(DF)
DT[, t.test(data=.SD, Age ~ Treated), by=Program]
   Program  statistic parameter   p.value   conf.int estimate null.value alternative
1: Program A -0.6286875  247.8390 0.5301326 -4.8110579 65.26667          0   two.sided
2: Program A -0.6286875  247.8390 0.5301326  2.4828527 66.43077          0   two.sided
3: Program B  1.4758524  230.5380 0.1413480 -0.9069634 67.15315          0   two.sided
4: Program B  1.4758524  230.5380 0.1413480  6.3211834 64.44604          0   two.sided
5: Program C  0.1994182  246.9302 0.8420998 -3.3560930 63.56557          0   two.sided
6: Program C  0.1994182  246.9302 0.8420998  4.1122406 63.18750          0   two.sided
7: Program D -1.1321569  246.0086 0.2586708 -6.1855837 62.31707          0   two.sided
8: Program D -1.1321569  246.0086 0.2586708  1.6701237 64.57480          0   two.sided
                method      data.name
1: Welch Two Sample t-test Age by Treated
2: Welch Two Sample t-test Age by Treated
3: Welch Two Sample t-test Age by Treated
4: Welch Two Sample t-test Age by Treated
5: Welch Two Sample t-test Age by Treated
6: Welch Two Sample t-test Age by Treated
7: Welch Two Sample t-test Age by Treated
8: Welch Two Sample t-test Age by Treated

In this format, for each Program, the statistic is the same for both and equal to t, the parameter here is the df, for conf.int, it goes (in order) lower then upper (so for Program A, the confidence interval is (-4.8110579, 2.4828527), and for estimate it will be group 0 and then group 1 (so for Program A, the mean for Treated == 0 is 65.26667, etc.

This was the quickest solution I could come up with, and you could loop through vars, or perhaps there's a simpler way.

---

EDIT: I only confirmed for Program A and for Age, using the following code:

DT[Program == 'Program A', t.test(Age ~ Treated)]
    Welch Two Sample t-test

data:  Age by Treated
t = -0.62869, df = 247.84, p-value = 0.5301
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.811058  2.482853
sample estimates:
mean in group 0 mean in group 1
       65.26667        66.43077

---

EDIT 2: Here is code that loops through your variables and rbind's them together:

do.call(rbind, lapply(vars, function(x) DT[, t.test(data=.SD, eval(parse(text=x)) ~ Treated), by=Program]))

Answer 2

You're getting errors because you're trying to access tt$statistic from within the function that creates tt. Some bracketing problems.

Here's one way to do it following your version

results <- lapply(progs, function (i) {
  DS = subset(DT, Program == i)
  o <- lapply(vars, function (i) {
    frm <- formula(paste0(i, '~ Treated'))
    tt <- t.test(frm, DS)
    data.frame(Variable=i, T=tt$statistic, P=tt$p.value)
  })
  o <- do.call(rbind, o)
  o$Program <- i
  o
})
do.call(rbind, results)

Or you can do it with rather rbind-ing using (e.g.) ddply (I think the rbinding still happens, just behind the scenes):

library(plyr)
combinations <- expand.grid(Program=progs, Y=vars)
ddply(combinations, .(Program, Y),
      function (x) {
          # x is a dataframe with the program and variable;
          #  just do the t-test and add the statistic & p-val to it
          frm <- formula(paste0(x$Y, '~ Treated'))
          tt <- t.test(frm, subset(DT, Program == x$Program))
          x$T <- tt$statistic
          x$P <- tt$p.value
          x
      })

Answer 3

You can get the same t-test out of a regression; if you think the effect of treatment is different for different programs, you should include an interaction. You can also specify multiple responses.

> m <- lm(cbind(Age,Dollars,DiseaseCnt)~Treated * Program - Treated - 1, DT)
> lapply(summary(m), `[[`, "coefficients")
$`Response Age`
                              Estimate  Std. Error       t value         Pr(>|t|)
ProgramProgram A         63.0875912409 1.294086510 48.7506752932 1.355786133e-265
ProgramProgram B         65.3846153846 1.400330869 46.6922616771 1.207761156e-252
ProgramProgram C         66.0695652174 1.412455172 46.7763979425 3.534894216e-253
ProgramProgram D         66.6691729323 1.313402302 50.7606640010 5.038015651e-278
Treated:ProgramProgram A  2.8593114140 1.924837595  1.4854819032  1.377339219e-01
Treated:ProgramProgram B -0.9786003470 1.919883369 -0.5097186438  6.103619649e-01
Treated:ProgramProgram C -0.5066022544 1.922108032 -0.2635659631  7.921691261e-01
Treated:ProgramProgram D -2.8657541289 1.919883369 -1.4926709484  1.358412980e-01

$`Response Dollars`
                                Estimate  Std. Error        t value     Pr(>|t|)
ProgramProgram A          998.5474452555 2.681598120 372.3702808887 0.0000000000
ProgramProgram B          997.4188034188 2.901757030 343.7292623810 0.0000000000
ProgramProgram C         1001.6869565217 2.926880936 342.2370019265 0.0000000000
ProgramProgram D         1001.2180451128 2.721624185 367.8752013053 0.0000000000
Treated:ProgramProgram A   -0.9899231316 3.988636646  -0.2481858388 0.8040419882
Treated:ProgramProgram B    2.5060086113 3.978370529   0.6299082986 0.5288996396
Treated:ProgramProgram C   -5.4721417069 3.982980462  -1.3738811324 0.1697889454
Treated:ProgramProgram D   -4.0043698991 3.978370529  -1.0065351806 0.3144036460

$`Response DiseaseCnt`
                               Estimate   Std. Error        t value         Pr(>|t|)
ProgramProgram A          4.53284671533 0.1793523653 25.27341475576 3.409326912e-109
ProgramProgram B          4.56410256410 0.1940771747 23.51694665775  1.515736580e-97
ProgramProgram C          4.25217391304 0.1957575279 21.72163675698  6.839384262e-86
ProgramProgram D          4.60150375940 0.1820294143 25.27890219412 3.133081901e-109
Treated:ProgramProgram A  0.13087009883 0.2667705543  0.49057175444  6.238378600e-01
Treated:ProgramProgram B -0.02274918064 0.2660839292 -0.08549625944  9.318841210e-01
Treated:ProgramProgram C  0.47375201288 0.2663922537  1.77840010867  7.564438017e-02
Treated:ProgramProgram D -0.31090546880 0.2660839292 -1.16844887901  2.429064705e-01

You specifically care about the Treated:Program entries of the regression table.