Five-ish Steps to Create Pretty Interaction Plots for a Multi-level Model in R

The present example uses intensive longitudinal data to examine how the effects of daily and average stressor exposure on negative affect may be buffered by daily and person-level control beliefs. We briefly run through preparatory steps and show the multi-level model used, then display how to plot the interaction effects. All in Five (ish) steps!

1. Preparation

Call in the necessary packages
- nlme, effects, and ggplot for running the model and making interaction plots; psych for descriptive statistics

library(nlme)
library(effects)

## Warning: package 'effects' was built under R version 3.3.1

library(ggplot2)
library(psych)

Read in the data
Check whether data are structured properly (e.g. binary, ordinal, or nominal variables are read as factors)

withscales <- read.csv("C:/Users/Rachel/Desktop/Berlin Summer/Stress and Control/withscales.csv")
scales <- read.csv("C:/Users/Rachel/Desktop/Berlin Summer/Stress and Control/scales.csv")
demog<-read.csv("C:/Users/Rachel/Desktop/Berlin Summer/Stress and Control/iSAHIB demographic variables_2015-05-27.csv")
#grand mean centering demographic data
summary(demog$age_may12_2010)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.98   30.16   48.18   47.10   57.59   88.71

summary(demog$edyears)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    2.00   15.00   16.50   16.36   18.00   24.00      18

Center variables and re-name variables, where necessary, to be more informative.

demog$c_age<-demog$age_may12_2010 - 47.10
demog$c_edu<-demog$edyears - 16.36
#re-labeling
withscales$DailyStressorOccurred<-withscales$stressorany_d
#making sure daily stressor and id are factors
withscales$DailyStressorOccurred<-as.factor(withscales$DailyStressorOccurred)

Merge data frames to create final data set

withscalesdemo<-Reduce(function(x, y) merge(x, y, by="id", all=TRUE), list(withscales,demog))
#make sure id is a factor variable
withscalesdemo$id<-as.factor(withscalesdemo$id)

2. Model time! Run the final model and look at the summary of the model to see what interactions you have.

For our example, this is the final model for which we want to plot the interactions.

lme15Sept<-lme(fixed= negaffect_d~
                   DailyStressorOccurred  +
                   PersMeanCent_DailyControl +  
                   DailyStressorOccurred*PersMeanCent_DailyControl +
                   
                   c_edu +
                   gender +
                   c_age +
                   GrandMeanCent_scale2 +  
                   GrandMeanCent_Exposure +  
                   GrandMeanCent_scale2*c_age + 
                   GrandMeanCent_Exposure*c_age + 
                   GrandMeanCent_Exposure*GrandMeanCent_scale2 +
                   
                   DailyStressorOccurred*c_age +
                   DailyStressorOccurred*GrandMeanCent_Exposure +
                   DailyStressorOccurred*GrandMeanCent_scale2 +
                   PersMeanCent_DailyControl*GrandMeanCent_scale2 +
                   PersMeanCent_DailyControl*c_age +
                   PersMeanCent_DailyControl*GrandMeanCent_Exposure +
                   
                   DailyStressorOccurred*PersMeanCent_DailyControl*GrandMeanCent_Exposure +
                   PersMeanCent_DailyControl*GrandMeanCent_scale2*c_age +
                   
                   DailyStressorOccurred*GrandMeanCent_scale2*c_age +
                   DailyStressorOccurred*PersMeanCent_DailyControl*c_age +
                   GrandMeanCent_Exposure*GrandMeanCent_scale2*c_age +
                   GrandMeanCent_Exposure*PersMeanCent_DailyControl*c_age,
                   
                   data=withscalesdemo, 
                   random = ~1 + DailyStressorOccurred |id,  
                   na.action= na.exclude)

Summary of the model! We won’t print this here, since it is extremely long for our model (shows all variable correlations).

summary(lme15Sept)

Instead, let’s look at the confidence intervals (CIs) for all fixed and random effects to find the interactions that do not contain 0 within the CI:

intervals(lme15Sept)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                                                                                 lower
## (Intercept)                                                              1.368088e+01
## DailyStressorOccurred1                                                   3.498318e+00
## PersMeanCent_DailyControl                                               -2.246312e-01
## c_edu                                                                   -4.158905e-01
## gender                                                                  -2.147970e+00
## c_age                                                                   -8.135979e-02
## GrandMeanCent_scale2                                                    -3.358713e+00
## GrandMeanCent_Exposure                                                   5.906586e+00
## DailyStressorOccurred1:PersMeanCent_DailyControl                        -1.348678e-01
## c_age:GrandMeanCent_scale2                                              -1.370424e-01
## c_age:GrandMeanCent_Exposure                                            -1.379742e-01
## GrandMeanCent_scale2:GrandMeanCent_Exposure                             -9.824631e+00
## DailyStressorOccurred1:c_age                                            -3.446900e-02
## DailyStressorOccurred1:GrandMeanCent_Exposure                           -3.165914e+00
## DailyStressorOccurred1:GrandMeanCent_scale2                             -1.735218e+00
## PersMeanCent_DailyControl:GrandMeanCent_scale2                          -1.132677e-04
## PersMeanCent_DailyControl:c_age                                         -6.585954e-04
## PersMeanCent_DailyControl:GrandMeanCent_Exposure                        -2.524181e-01
## DailyStressorOccurred1:PersMeanCent_DailyControl:GrandMeanCent_Exposure -3.651457e-01
## PersMeanCent_DailyControl:c_age:GrandMeanCent_scale2                    -2.220832e-03
## DailyStressorOccurred1:c_age:GrandMeanCent_scale2                       -5.166073e-02
## DailyStressorOccurred1:PersMeanCent_DailyControl:c_age                  -5.068820e-03
## c_age:GrandMeanCent_scale2:GrandMeanCent_Exposure                       -3.456495e-01
## PersMeanCent_DailyControl:c_age:GrandMeanCent_Exposure                  -3.668922e-05
##                                                                                 est.
## (Intercept)                                                             15.483998899
## DailyStressorOccurred1                                                   4.328499731
## PersMeanCent_DailyControl                                               -0.179139413
## c_edu                                                                   -0.105096954
## gender                                                                   0.292982588
## c_age                                                                   -0.006197596
## GrandMeanCent_scale2                                                    -1.996370276
## GrandMeanCent_Exposure                                                  11.644922779
## DailyStressorOccurred1:PersMeanCent_DailyControl                        -0.085902780
## c_age:GrandMeanCent_scale2                                              -0.060289902
## c_age:GrandMeanCent_Exposure                                             0.103516970
## GrandMeanCent_scale2:GrandMeanCent_Exposure                             -4.368763898
## DailyStressorOccurred1:c_age                                             0.009259742
## DailyStressorOccurred1:GrandMeanCent_Exposure                            0.091880064
## DailyStressorOccurred1:GrandMeanCent_scale2                             -0.897350560
## PersMeanCent_DailyControl:GrandMeanCent_scale2                           0.014348318
## PersMeanCent_DailyControl:c_age                                          0.001466396
## PersMeanCent_DailyControl:GrandMeanCent_Exposure                        -0.109715573
## DailyStressorOccurred1:PersMeanCent_DailyControl:GrandMeanCent_Exposure -0.194117580
## PersMeanCent_DailyControl:c_age:GrandMeanCent_scale2                    -0.001387604
## DailyStressorOccurred1:c_age:GrandMeanCent_scale2                       -0.002135321
## DailyStressorOccurred1:PersMeanCent_DailyControl:c_age                  -0.002635248
## c_age:GrandMeanCent_scale2:GrandMeanCent_Exposure                       -0.086312985
## PersMeanCent_DailyControl:c_age:GrandMeanCent_Exposure                   0.003719078
##                                                                                 upper
## (Intercept)                                                             17.2871174553
## DailyStressorOccurred1                                                   5.1586818152
## PersMeanCent_DailyControl                                               -0.1336476308
## c_edu                                                                    0.2056965807
## gender                                                                   2.7339348321
## c_age                                                                    0.0689645952
## GrandMeanCent_scale2                                                    -0.6340276049
## GrandMeanCent_Exposure                                                  17.3832598131
## DailyStressorOccurred1:PersMeanCent_DailyControl                        -0.0369378022
## c_age:GrandMeanCent_scale2                                               0.0164626060
## c_age:GrandMeanCent_Exposure                                             0.3450081672
## GrandMeanCent_scale2:GrandMeanCent_Exposure                              1.0871036844
## DailyStressorOccurred1:c_age                                             0.0529884878
## DailyStressorOccurred1:GrandMeanCent_Exposure                            3.3496745977
## DailyStressorOccurred1:GrandMeanCent_scale2                             -0.0594829013
## PersMeanCent_DailyControl:GrandMeanCent_scale2                           0.0288099041
## PersMeanCent_DailyControl:c_age                                          0.0035913871
## PersMeanCent_DailyControl:GrandMeanCent_Exposure                         0.0329869875
## DailyStressorOccurred1:PersMeanCent_DailyControl:GrandMeanCent_Exposure -0.0230894667
## PersMeanCent_DailyControl:c_age:GrandMeanCent_scale2                    -0.0005543766
## DailyStressorOccurred1:c_age:GrandMeanCent_scale2                        0.0473900863
## DailyStressorOccurred1:PersMeanCent_DailyControl:c_age                  -0.0002016760
## c_age:GrandMeanCent_scale2:GrandMeanCent_Exposure                        0.1730235661
## PersMeanCent_DailyControl:c_age:GrandMeanCent_Exposure                   0.0074748447
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                                              lower        est.     upper
## sd((Intercept))                          5.4182562  6.47752016 7.7438692
## sd(DailyStressorOccurred1)               1.9239372  2.68887896 3.7579553
## cor((Intercept),DailyStressorOccurred1) -0.5404973 -0.06000997 0.4499946
## 
##  Within-group standard error:
##    lower     est.    upper 
## 8.833355 8.973700 9.116276

We find that there are two 2-way interactions and three three 3-way interactions for us to plot!

3. Obtain fitted values for the interaction of interest

Beginning with the first 2-way interaction effect (DailyStressorOccurred*PersMeanCent_DailyControl), use the effects() function from the effects package to give us information about the fitted values, standard errors, and lower and upper confidence intervals for the values of the interaction variables.
- In the argument term=, specify the interaction (exactly as it is written in the model or results of the model)
- If you don’t want the function to automatically select the values of your predictors, then use the the argument xlevels= to specify vectors of values you would like to set for each numeric predictor. If you would like to use one standard deviation (SD) above and below the mean, use the describe function for the variable of interest from your original data frame. Then use those values in the xlevels= argument as a list of vectors: xlevels= list(predictor1=c(-1SD, 1SD), predictor2=c(-1SD, 1SD)).
  - Note: For factor predictors, the effects function will include all levels of the factor. You may only specify factor levels by changing the factor predictor to a numeric predictor prior to running the model in nlme.)
- In the argument mod=, specify the variable name for the model from which you are taking this interaction
- Put the effects object into a variable, so you can convert it to a data frame and use it in your plots!

describe(withscalesdemo$PersMeanCent_DailyControl) #for standard deviations

##   vars    n mean    sd median trimmed  mad   min   max  range  skew
## 1    1 8533    0 13.71    1.4    0.91 9.86 -95.6 58.82 154.43 -1.15
##   kurtosis   se
## 1     4.67 0.15

ef1 <- effect(term="DailyStressorOccurred*PersMeanCent_DailyControl", xlevels= list(PersMeanCent_DailyControl=c(-13.71, 13.71)), mod=lme15Sept)

## NOTE: DailyStressorOccurred:PersMeanCent_DailyControl is not a high-order term in the model

efdata1<-as.data.frame(ef1) #convert the effects list to a data frame
efdata1 #print effects data frame

##   DailyStressorOccurred PersMeanCent_DailyControl      fit        se
## 1                     0                    -13.71 17.80431 0.7947768
## 2                     1                    -13.71 23.28865 0.6951044
## 3                     0                     13.71 12.98323 0.7405723
## 4                     1                     13.71 16.10892 0.7023356
##      lower    upper
## 1 16.24634 19.36228
## 2 21.92606 24.65124
## 3 11.53151 14.43494
## 4 14.73216 17.48568

Note, the effects() function will warn you when the interaction you are looking at is superceded by a higher order interaction. This does not mean you should not plot it. It just means that you need to be aware of the higher order interactions which form the context for the lower-order one you are plotting.

4. Prepare effects data to make appropriate plot

We will recode Daily Control Beliefs (PersMeanCent_DailyControl) as a factor variable to make it easier to use in our plot.

efdata1$PersMeanCent_DailyControl<-as.factor(efdata1$PersMeanCent_DailyControl)

5. Now to make the actual interaction plots!

5a) Two-way interaction

DailyStressorOccurred*PersMeanCent_DailyControl

Use ggplot() to create a plot layer. Within this layer, the first argument is the data frame you are using (the one you created from the effects list), and the second argument is the aesthetics aes(), in which you specify the following:
- x= is where you write in your main independent variable of interest (let’s call this IV1), which will be plotted across the x-axis
- y= is where you write fit (which was automatically generated and included in the effects data frame), so you will plot the fitted values of IV1 across the y-axis
- color= is where you put the independent variable (IV2) that moderates the effects of IV1. This produces groupings of IV1’s fitted values for each unique value of IV2
- group= is an additional option where you can insert IV2 if you would like error bands instead of error bars that we add in another layer.
Add a geom_point() layer to plot the above x and y values as points
Add a geom_line() layer to draw a line connecting the above fit points for each group. Note: If you wanted error bars (which we add in the geom_ribbon layer) instead of bands, add an aesthetic here that specifies group=PersonMeanCent_DailyControl.
Add a geom_ribbon layer to add error bars/bands. You need to specify the min and max of the error bars, by using the standard error variable (se) from the effects data frame. So, in aes() here, you specify ymin= fit - se, ymax= fit + se. *Then you’re all set to label your axes and give a plot title, and you’re off!

ggplot(efdata1, aes(x=DailyStressorOccurred, y=fit, color=PersMeanCent_DailyControl,group=PersMeanCent_DailyControl)) + 
    geom_point() + 
    geom_line(size=1.2) +
    geom_ribbon(aes(ymin=fit-se, ymax=fit+se, fill=PersMeanCent_DailyControl),alpha=0.3) + 
    labs(title = "Reactivity by Daily Control Beliefs", x= "Daily Stressor (0=Did not occur, 1=Occurred)", y="Negative Affect", color="Daily Control Beliefs", fill="Daily Control Beliefs") + theme_classic() + theme(text=element_text(size=20))

5b) Three-way interaction:

DailyStressorOccurred1:PersMeanCent_DailyControl:GrandMeanCent_Exposure

First, we’ll follow steps 3 and 4 for this interaction effect.

Step 3 Obtain fitted values for the interaction of interest

describe(withscalesdemo$PersMeanCent_DailyControl)

##   vars    n mean    sd median trimmed  mad   min   max  range  skew
## 1    1 8533    0 13.71    1.4    0.91 9.86 -95.6 58.82 154.43 -1.15
##   kurtosis   se
## 1     4.67 0.15

describe(withscalesdemo$GrandMeanCent_Exposure)

##   vars    n  mean   sd median trimmed  mad   min  max range  skew kurtosis
## 1    1 8558 -0.01 0.28   0.08    0.03 0.24 -0.72 0.28     1 -0.94    -0.25
##   se
## 1  0

ef3 <- effect(term="DailyStressorOccurred*PersMeanCent_DailyControl*GrandMeanCent_Exposure", xlevels=list(PersMeanCent_DailyControl=c(-13.71, 13.71), GrandMeanCent_Exposure= c(-.28, .28)) , mod=lme15Sept)
efdata3<-as.data.frame(ef3)
summary(efdata3)

##  DailyStressorOccurred PersMeanCent_DailyControl GrandMeanCent_Exposure
##  0 :4                  Min.   :-13.71            Min.   :-0.28         
##  1 :4                  1st Qu.:-13.71            1st Qu.:-0.28         
##  NA:0                  Median :  0.00            Median : 0.00         
##                        Mean   :  0.00            Mean   : 0.00         
##                        3rd Qu.: 13.71            3rd Qu.: 0.28         
##                        Max.   : 13.71            Max.   : 0.28         
##       fit              se             lower            upper      
##  Min.   :10.28   Min.   :0.9862   Min.   : 8.174   Min.   :12.38  
##  1st Qu.:14.26   1st Qu.:1.0518   1st Qu.:12.091   1st Qu.:16.44  
##  Median :17.14   Median :1.1130   Median :15.053   Median :19.23  
##  Mean   :17.70   Mean   :1.1015   Mean   :15.539   Mean   :19.86  
##  3rd Qu.:19.71   3rd Qu.:1.1451   3rd Qu.:17.435   3rd Qu.:21.99  
##  Max.   :27.92   Max.   :1.2387   Max.   :25.985   Max.   :29.85

Step 4 Prepare effects data to make appropriate plot

This step is slightly different for a three-way interaction than for a two-way, as we now must incorporate a third independent variable. We still want to make IV2 into a factor type variable.
But now we also want to create a factor for IV3, and break IV3 into the groups by which we want to separate our plots (i.e. low vs high).
- To aid our thinking, we first create a variable facets which gets a vector of the labels “low” and “high”.
- Then we use the cut() function, which will allow us to cut IV3’s numeric data into bins for a factor variable. Within this function, we put the arguments:
  - data= the data we are using (in this case, efdata3$GrandMeanCent_Exposure)
  - breaks= the breaks, which are the vector of cutpoints used to partition the numeric data into bins (in this case, we have two bins, defined by the breaks. The data is cut from -0.2 to 0 for one bin, and 0 to 0.2 for the second bin).
  - labels= the labels for the bins (which we conveniently had created earlier as our facets variable)
  - ordered_result= an indication of whether this re-defined factor variable is ordinal (i.e. does it matter than low is less than high. In our case, yes, to be appropriate in our data conceptualization.)
We then use the subset command to re-code our data frame, allowing us to use only the rows in which our IV3 contains the facet values (“low” or “high”)

efdata3$PersMeanCent_DailyControl<-factor(efdata3$PersMeanCent_DailyControl) #making IV2 a factor, as usual
efdata3$DailyStressorOccurred<-factor(efdata3$DailyStressorOccurred)

Step 5 Now to make the actual interaction plots!

Use ggplot() to create a plot layer. Within this layer, the first argument is the data frame you are using (the one you just subsetted from your original effects data frame), and the second argument is the aesthetics aes(), in which you specify the following:
- x= is where you write in your main independent variable of interest (let’s call this IV1), which will be plotted across the x-axis
- y= is where you write fit (which was automatically generated and included in the effects data frame), so you will plot the fitted values of IV1 across the y-axis
- color= is where you put the independent variable (IV2) that moderates the effects of IV1. This produces groupings of IV1’s fitted values for each unique value of IV2
- In this case, we will add an aesthetic here that specifies group=PersonMeanCent_DailyControl, so our error bands (which we add in the geom_ribbon layer) will be separated by group.
Add a geom_point() layer to plot the above as points
Add a geom_line() layer to draw a line connecting the above fit points for each group.
Add a geom_ribbon layer to add error bars/bands. You need to specify the min and max of the error bars, by using the standard error variable (se) from the effects data frame. So, in aes() here, you specify ymin= fit - se, ymax= fit + se.
The only new step here as opposed to our 2-way interaction is that we now want to create two plots that show us the above 2-way interaction in the presence of low exposure or high exposure separately. We do this by adding a facet_wrap layer in which we specify the groupping variable by putting a ~ in front of it (e.g. facet_wrap(~GrandMeanCent_Exposure)) *Then you’re all set to label your axes and give a plot title, and you’re off! (Note the two separate plots are already labeled the way you specified in your facets variable)

ggplot(data=efdata3, aes(x=DailyStressorOccurred, y=fit, group=PersMeanCent_DailyControl)) + geom_point(aes(color=PersMeanCent_DailyControl)) + geom_line(aes(color=PersMeanCent_DailyControl))  + 
    geom_ribbon(aes(ymin=fit-se,ymax=fit+se, fill=PersMeanCent_DailyControl),alpha=0.3) + facet_wrap(~GrandMeanCent_Exposure) + labs(title = "Stressor Reactivity by \n Daily Control Beliefs \n Facetted by Exposure Groups", x= "Daily Stressor (0=Did not occur, 1=Occurred)", y="Negative Affect", color="Daily Control Beliefs", fill="Daily Control Beliefs")  + theme_classic() + theme(text=element_text(size=20))

Additional Examples

DailyStressorOccurred1:GrandMeanCent_scale2

describe(withscalesdemo$GrandMeanCent_scale2)

##   vars    n mean  sd median trimmed  mad   min  max range  skew kurtosis
## 1    1 8558 0.02 1.1   0.05    0.06 1.24 -2.89 1.83  4.72 -0.35     -0.4
##     se
## 1 0.01

ef4 <- effect(term="DailyStressorOccurred:GrandMeanCent_scale2", xlevels=list(GrandMeanCent_scale2=c(-1.1,1.1)), mod=lme15Sept)

## NOTE: DailyStressorOccurred:GrandMeanCent_scale2 is not a high-order term in the model

efdata4<-as.data.frame(ef4)


efdata4$DailyStressorOccurred<-as.factor(efdata4$DailyStressorOccurred)
efdata4$GrandMeanCent_scale2<-as.factor(efdata4$GrandMeanCent_scale2)

ggplot(efdata4, aes(x=DailyStressorOccurred, y=fit, color=GrandMeanCent_scale2, group=GrandMeanCent_scale2)) + geom_point() + geom_line(size=1.2) +
    geom_ribbon(aes(ymin=fit-se,ymax=fit+se, fill=GrandMeanCent_scale2),alpha=0.3) + labs(title = "Daily Stressor Occurrence by General Control Beliefs", x= "Daily Stressor Occurrence", y="Negative Affect", color="Personal Control Beliefs", fill="Personal Control Beliefs") + theme_classic() + theme(text=element_text(size=20))

PersMeanCent_DailyControl:c_age:GrandMeanCent_scale2

describe(withscalesdemo$c_age)

##   vars    n mean    sd median trimmed   mad    min   max range skew
## 1    1 8558 0.84 18.52   2.27    0.06 23.29 -28.12 41.61 69.73 0.25
##   kurtosis  se
## 1    -0.81 0.2

describe(withscalesdemo$PersMeanCent_DailyControl)

##   vars    n mean    sd median trimmed  mad   min   max  range  skew
## 1    1 8533    0 13.71    1.4    0.91 9.86 -95.6 58.82 154.43 -1.15
##   kurtosis   se
## 1     4.67 0.15

describe(withscalesdemo$GrandMeanCent_scale2)

##   vars    n mean  sd median trimmed  mad   min  max range  skew kurtosis
## 1    1 8558 0.02 1.1   0.05    0.06 1.24 -2.89 1.83  4.72 -0.35     -0.4
##     se
## 1 0.01

ef5 <- effect(term="PersMeanCent_DailyControl:c_age:GrandMeanCent_scale2",xlevels=list(c_age=c(-18.52, 18.52),PersMeanCent_DailyControl=c(-13.71,13.71), GrandMeanCent_scale2=c(-1.1,1.1)), mod=lme15Sept)
efdata5<-as.data.frame(ef5)

efdata5$PersMeanCent_DailyControl<-as.factor(efdata5$PersMeanCent_DailyControl)
efdata5$GrandMeanCent_scale2<-as.factor(efdata5$GrandMeanCent_scale2)
efdata5$c_age<-as.factor(efdata5$c_age)

ggplot(efdata5, aes(x=PersMeanCent_DailyControl, y=fit, color=GrandMeanCent_scale2, group=GrandMeanCent_scale2, fill=GrandMeanCent_scale2)) + geom_point() + geom_line(size=1.2) +
    geom_ribbon(aes(ymin=fit-se,ymax=fit+se),alpha=0.3) + facet_wrap( ~c_age) + labs(title = "Daily Control Beliefs by General Control Beliefs Facetted by Age", x= "Daily Control Beliefs", y="Negative Affect", color="Personal Control Beliefs", fill="Personal Control Beliefs") + theme_classic() +theme(text=element_text(size=20))

DailyStressorOccurred1:PersMeanCent_DailyControl:c_age

ef6 <- effect(term="DailyStressorOccurred:PersMeanCent_DailyControl:c_age",xlevels=list(PersMeanCent_DailyControl=c(-13.71,13.71),c_age=c(-18.52,18.52)), mod=lme15Sept)
efdata6<-as.data.frame(ef6)
efdata6$PersMeanCent_DailyControl<-as.factor(efdata6$PersMeanCent_DailyControl)
efdata6$c_age<-as.factor(efdata6$c_age)
ggplot(efdata6, aes(x=DailyStressorOccurred, y=fit, color=PersMeanCent_DailyControl, group=PersMeanCent_DailyControl, fill=PersMeanCent_DailyControl)) + geom_point() + geom_line(size=1.2) +
    geom_ribbon(aes(ymin=fit-se,ymax=fit+se),alpha=0.3) + facet_wrap( ~c_age) + labs(title = "Stressor Reactivity by Daily Control Beliefs Facetted by Age", x= "Stressor Occurrence (0=Did not occur, 1=Occurred)", y="Negative Affect", color="Daily Control Beliefs", fill="Daily Control Beliefs") + theme_classic() +theme(text=element_text(size=20))