Revisiting the t-test: just another linear model

Before I even start with this topic, I want to point out that I took much of the information I present here from Jonas Kristoffer Lindelov’s fantastic blog post “Common statistical tests are linear models (or: how to teach stats)”. If you have not read it yet, you should definitely do so, either now or after you finished this part of the tutorial. Seriously, it is a great blog-post that inspired me to write up this part of the tutorial in a similar way in which we want to learn one flexible method of data-simulation rather than changing what we are doing for each different situation.

Ok so what does it mean to say that, for instance, a t-test is just a linear model? First of all, here is a quick reminder about the linear model with a quick example:

\[ y_i = \beta_0 + \beta_1x_i \]

The above formula states that each value of y (the little i can be read as “each”) can be described by adding a value \(\beta_0\) and the product of \(beta_0 \times x\). Now, how does this relate to the t-test? Again, if you want a more in-depth discussion read the blog-post above. To demonstrate the point, lets revisit the example from part II of this tutorial. We had 2 groups and we assumed that one group can be described as group1 = normal(0,2) and the other group can be described as group2 = normal(1,2). Now if we want to test whether the difference in means is significant, we can run a independent-samples t-test. However, if we rephrase the problem what we want to know whether group membership (whether you belong to group1 or group2; the \(x\) variable in this case) significantly predicts the person’s score (y in this case). To do so, we can just assume let \(beta_0\) describe the mean of, say, group1 in this case. If we do this and if we re-write model above to predict the mean of group1, it looks the following way:

\[ M_{group1} = M_{group1}+\beta_1\times x \]

What should \(\beta_1\) and \(x\) be in this case? Well, as we already have the mean of group1 on the right side in the term that was previously \(beta_0\), we do not want to add anything to it. Thus, if someone member from group1, to predict the mean-score we should have an observed value for group-membership of 0. This way, the formula becomes

\[ M_{group1} = M_{group1}+\beta_1\times 0 \]

and the \(\beta_1\) term becomes 0 entirely. Note that we of course do not know the mean of group1 as we will simulate it. \(\beta_0\) (or \(M_{group1}\) here) is a parameter, i.e. something that we want the model to estimate because we do not know it.

What if we want to describe the mean of a person in group2? We can do this by just adding something to the mean of group1, namely the difference between the two groups. In other words, we can describe the mean of group2 by saying its the mean of group1 + the difference of their means times x:

\[ M_{group2} = M_{group1}+(M_{group2}-M_{group1})\times x \]

What should \(\beta_1\) and \(x\) be now? Just like \(\beta_0\), our new mean-difference parameter \(\beta_1\) is an unknown that depends on our simulated means and will be estimated by the model. \(x\) on the other hand should be 1, because if someone is from group2 we want the predicted mean of that person to be the mean of group1 + 1 times the difference between group1 and group2.

Therefore if someone is from group2 our model becomes:

\[ M_{group2} = M_{group1}+(M_{group2}-M_{group1})\times 1 \]

What we just did is so-called dummy coding, i.e. we made a new dummy variable that tells us someone’s group membership:

• 0 if the person is in group1
• 1 if the person is in group2

If we use this dummy variable in the linear model we can do exactly the same test that we did with the t-test earlier. Let me demonstrate with the exact same groups that we simulated in the beginning of part II of this tutorial.

set.seed(1234)
group1 <- rnorm(30, 1, 2)
group2 <- rnorm(30, 0, 2)

Now we can run the independent-samples t-test again, yielding the same results as in part II.

t.test(group1, group2, paired = FALSE, var.equal = TRUE, conf.level = 0.9)

## 
##  Two Sample t-test
## 
## data:  group1 and group2
## t = 3.1402, df = 58, p-value = 0.002656
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
##  0.7064042 2.3143690
## sample estimates:
## mean of x mean of y 
##  0.407150 -1.103237

To do our linear model analysis instead, we first have to create a data-frame that contains the 2 groups and tells us about their group membership. In other words we will make a data-set that has the variables score which are the values in the vectors of group1 and group2 and a variable group that will be 0 if a person is from group1 and 1 if the person is from group2. In this case, after putting the data together, we create a dummy variable with ifelse a function that you can read like: if the group variable has value “group1” in a certain row, then make dummy_group in that row = 0, else make it 1.

names(group1) <- rep("group1", length(group1)) # name the vectors to use the names as the group variable
names(group2) <- rep("group2", length(group2))

lm_groupdat <- data.frame(score = c(group1,group2), group = c(names(group1), names(group2))) # make data-set from scores and names

lm_groupdat$dummy_group <- ifelse(lm_groupdat$group == "group1", 0, 1) # create dummy variable 

Now we can use this to run our linear model with the lm function. In case you have not used it before, we indicate the intercept, i.e. \(beta_0\) with a 1 in this case, to tell the model that we want 1 times \(beta_0\) in there and tell it that we want dummy_group times the difference between the groups, \(beta_1\) in there. The code looks like this:

summary(lm(score ~ 1+ dummy_group, dat = lm_groupdat)) # use summary function to get p-values

## 
## Call:
## lm(formula = score ~ 1 + dummy_group, data = lm_groupdat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0985 -1.1327 -0.4088  0.7360  5.4245 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.4072     0.3401   1.197  0.23612   
## dummy_group  -1.5104     0.4810  -3.140  0.00266 **
## 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

We see that, if we collect 10,000 people per group, everything would be significant with focus having the biggest effect (look at the ges column in the table for generalized eta-squared values) and with media having the smallest effect (in fact close to 0 as we wanted). That even the main-effect for media is still significant just shows that this happens even with very small deviations from the exact null-hypothesis in these huge samples.

Note that again we could rewrite this anova as a dummy-coded linear model as:

\[ score_i = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 \]

By making dummy-variables we can code this in lm

aov_dat$focus_dummy <- ifelse(aov_dat$focus == "internal", 0 , 1)
aov_dat$media_dummy <- ifelse(aov_dat$media == "text", 0 , 1)

lm_int <- lm(score ~ 1 + focus_dummy + media_dummy + focus_dummy:media_dummy, data = aov_dat)

From the estimates that we get here we see how the model is estimated:

\(\beta_0\) = mean(group_TI) as both dummies are 0 for this one. \(\beta_1\) = mean(group_TE)-mean(group_TI) as the focus variable changes from 0 to 1 within the text-group \(\beta_2\) = mean(group_VI)-mean(group_TI) as the media variable changes from 0 to 1 within the internal focus group \(\beta_3\) = mean(group_VE)-mean(group_TE)+mean(group_TI)-mean(group_VI) to see what a crossover of the variables (i.e. the interaction) compares to \(\beta_0\)

If we want to test whether the effect of the interaction, for instance is significant, we can conduct an ANOVA for the interaction to see how the model with the interaction compares to a model without the interaction, which we do in the following code-block by fitting a model that does not have the interaction in their and compare it with the anova function.

lm_null <- lm(score ~ 1 + focus_dummy + media_dummy, data = aov_dat)
anova(lm_null, lm_int)

## Analysis of Variance Table
## 
## Model 1: score ~ 1 + focus_dummy + media_dummy
## Model 2: score ~ 1 + focus_dummy + media_dummy + focus_dummy:media_dummy
##   Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
## 1 399997 90721027                                  
## 2 399996 90112027  1    609000 2703.3 < 2.2e-16 ***
## 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We see immediately that the estimates changed somewhat. Now, for example the intercept represents the mean of all 4 other means (i.e. the grand mean):

mean(c(mean(group_TI), mean(group_TE), mean(group_VE), mean(group_VI))) = 43.775609.

and the first estimate represents the difference of each focus level from that grand-mean, for the “internal groups this is”

mean(c(mean(group_TI), mean(group_VI)))-mean(c(mean(group_TI), mean(group_TE), mean(group_VE), mean(group_VI))) = -2.2276811

and for the external groups this is

mean(c(mean(group_TE), mean(group_VE)))-mean(c(mean(group_TI), mean(group_TE), mean(group_VE), mean(group_VI))) = 2.2276811

respectively.

Why would we use this coding-scheme in the linear model instead of the other one? The p-value for the interaction is exactly the same in both cases. However, if we compare the p-values of main-effects here, we will see how they differ.

For instance, if we tell aov_car that we want to keep our dummy-coding instead of changing it, we can see how the two outputs differ.

aov_car(score ~ focus*media+ Error(participant), data = aov_dat, type = 3, afex_options("check_contrasts" = FALSE))

## Anova Table (Type 3 tests)
## 
## Response: score
##        Effect        df    MSE           F  ges p.value
## 1       focus 1, 399996 225.28  876.77 *** .002  <.0001
## 2       media 1, 399996 225.28  829.66 *** .002  <.0001
## 3 focus:media 1, 399996 225.28 2703.27 *** .007  <.0001
## 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.94 on 3992 degrees of freedom
## Multiple R-squared:  0.03674,    Adjusted R-squared:  0.03505 
## F-statistic: 21.75 on 7 and 3992 DF,  p-value: < 2.2e-16

Note how the effects that we put in before are still there but they are weaker now as we have smokers included in our sample. Moreover, smoking causes a more positive attitude towards smoking, just as we wanted, and the interaction between media and focus is not influenced by whether someone smokes or not, as we put the same punishment for the treatment effect in all the groups, leaving the interaction uninfluenced by whether someone smokes or not.

Summary

So now we have seen how we can simulate data for analyses that many people use in their everyday life as a (psychological) researcher. So far, all of the power-analyses that we did can also be done in standard software like G*Power (though in part I of this tutorial I explained why I still prefer simulation). However, as nice as these analyses are, we often deal with more complicated situations where we have to deal with hierarchical or cross-classified data. In these cases, other approaches like linear-mixed effects models, also called hierarchical models, need to be used to analyse our data. This is where simulation really shines as the structure of the data strongly depends on the research design and in a point-and-click software (though existing) it is difficult to offer all the flexibility that we might need. Thus, in the fourth part, we will see how we can simulate data and run a power-analysis in these cases.