Simulation as a tool for personality research

Prepared as an Appendix to Experimental Approaches to the Study of Personality
William Revelle
Northwestern University

To appear in: Personality Research Methods, edited by Robins, Fraley and Krueger
Contents:
  1. Using the R computer package as a tool for personality research
  2. Data simulation as a research tool
  3. Data organization and simple summary statistics
  4. Using the general linear model to fit data
  5. Graphical presentations techniques
  6. Scaling artifacts -- some graphical explorations

Using the R computer package as a tool for personality research

It is useful when learning how to do personality research to consider what ideal results would look like and also to consider the problem of various statistical artifacts. This page, to accompany a chapter on Experimental Approaches to Personality, provides R code for several different simulated studies. For an introduction on how to use R for personality research, consult my notes for personality research. For a somewhat more extensive introduction, try my tutorial for R for Psychological Research.

There are many possible statistical programs that can be used in psychological research. They differ in multiple ways, at least some of which are ease of use, generality, and cost. Some of the more common packages used are Systat, SPSS, and SAS. These programs have GUIs (Graphical User Interfaces) that are relatively easy to use but that are unique to each package. These programs are also very expensive and limited in what they can do. Although convenient to use, GUI based operations are difficult to discuss in written form. When teaching statistics or communicating results to others, it is helpful to use examples that may be used by others using different computing environments and perhaps using different software. The simulations and analyses presented here use an alternative approach that is widely used by practicing statisticians, the statistical environment R. This is not meant as a user's guide to R, but merely a demonstration of the power of R to simulate and analyse personality data. I hope that enough information is provided in this brief guide to make the reader want to learn more. (For the impatient, a brief guide to analyzing data for personality research is also available.)

It has been claimed that "The statistical programming language and computing environment S has become the de-facto standard among statisticians. The S language has two major implementations: the commercial product S-PLUS, and the free, open-source R. Both are available for Windows and Unix/Linux systems; R, in addition, runs on Macintoshes." (From John Fox's short course on S and R).

Data simulation as a research tool

Simulation of data is a surprisingly powerful research tool. For by forcing the theorist to specify the theoretical predictions precisely enough to generate artificial data, theoretical issues are clarified. Parameters need to be specified and relationships need to be identified, not by mere hand waving, but rather by accurate identification. The amount of error in measurement needs to be directly specified in the simulation, requiring the researcher to think about issues of reliability and the shape of latent-observed relationships.

In addition to forcing the theorist to specify the model, by knowing what the right answer is, one can test the power of one's analytical tools to detect the correct model. By knowing the true answer, and observing that small samples do not necessarily detect effects that are there (Type II errors) and that by doing many analyses, the probability of detecting effects that are not there (Type I errors) is higher than expected, the researcher learns greater respect for the problems of good design and good analysis.

By showing the R code for generating simulated data as well as the code for analysis the reader is encouraged to try out some of the techniques for themselves rather than just passively reading yet another chapter.

In all the simulations, data are generated based upon a model that Positive Affect is an interactive and non-linear function of Extraversion and rewarding cues in the environment, that Negative Affect is an additive but non-linear function of Neuroticism and punishing cues in the environment, that arousal is an additive function of stimulant drugs and Introversion, and that cognitive performance is a curvilinear function of arousal. This model is the “truth” that the analyses hope to uncover. Unlike real studies, in all simulations we have access to the latent (but unobserved) “true” variables as well as the observed Person and Outcome variables. (The simulations are based upon somewhat simplified summaries of a number of studies conducted over the past years at the Personality, Motivation and Cognition lab at Northwestern University but are idealizations and simplifications of the actual results of those studies.) Readers are encouraged to try the simulations and analyses for themselves, varying the sample sizes and the strength of the relationships. By presetting the seed for the random number generator to a memorable constant value (Adams, 1980) the results obtained in these simulations and reported below should match those carried out by copying the code in the appendix and running the simulations.

The first section is a simulation of how the personality dimensions of Introversion-Extraversion and Emotional Stability- Neuroticism interact with rewards and punishments to affect Positive and Negative Affect. The simulation uses particular parameter settings that were chosen for pedagogical purposes and only partly reflect published results. For a more detailed discussion of the psychological models being simulated, consult the chapter. For an introduction to personality theory and research, visit the Personality-Project.

The second simulation models the interactive effect on cognitive performance of the personality variable introversion-extraversion with arousal inductions such as caffeine. This is also based upon work at the Personality, Motivation and Cognition lab.

The third simulation is concerned with problems of scaling and shows how it important it is to choose the appropriate metric.

All simulations make use of some useful additions to R that may be found in the "psych" package at the personality-project.org. To use these additions, first load R and then use the package manager and installer.

Simulation 1: Personality and Rewards and Punishments

The following code is directly executable in R. The user is encouraged to modify the parameters to explore further how to detect main effects and interactions in personality research. The text in blue is directly executable and should be copied into R. Results and graphics from these commands are included for comparison purposes.

Note: To reproduce the same "random" simulations each time, include the set.seed(42) command. To make a different simulation each time, do not include that line.


#R code for simulating personality x situation effects
#
#Part 1: Simulate the data

set.seed(42)   #random number seed is fixed to produce identical "random" sequences
               #remove to allow each run to be different
               
#first set some parameters of the model
#change these to try different models

num <- 100    		#number of people to simulate
weightreward <- .5 #an index of how strong is the effect of reward on positive affect
weightpunish <- .4  #how strong is the effect of punishment on negative affect
weight_e<-  .0 		#how strong is the effect of extraversion on positive affect
weight_n<-  .3 		#how strong is the effect of neuroticism on negative affect?
weight_er <-.4 		#how strong is the interaction of e x reward on positive affect
weight_np <- 0 		#how strong is the interactionof n * punish on negative affect
reliability_e <- .7	#the reliability of the observed extraversion score
reliability_n <-.8   #the reliability of the observed neuroticism score
reliability_P <-.7   #the reliability of observed measures of Postive Affect
reliability_N <- .8  #the reliability of observed measures of Negative Affect
weight_arousal <- .7	#relationship between extraversion and arousal
reliability_arousal <- .8 #reliability of arousal
weight_caff <- .5    #within subject weight of effect of caffeine
mean_E <- 3 		#mean of true extraversion
mean_N <- 3			#mean of true neuroticism

#generate the data using the random normal distribution 

						#first simulate true (latent) scores
true_e <- rnorm(num,mean_E)   #true trait extraversion is normally distributed with sigma=1
true_n <- rnorm(num,mean_N)   	#true trait neuroticism is normally distributed 
true_arousal_plac <- rnorm(num) - weight_arousal * (true_e - mean_E) -weight_caff
true_arousal_caff <- rnorm(num) - weight_arousal * (true_e - mean_E) + weight_caff
                        #observed E and N are a mixture of true and error scores
extraversion <- sqrt(reliability_e)*(true_e-mean_E) + sqrt(1-reliability_e)*rnorm(num) + mean_E
neuroticism <- sqrt(reliability_n)*(true_n -mean_N)+ sqrt(1-reliability_n)*rnorm(num) + mean_N
arousal_plac <- sqrt(reliability_arousal) * (true_arousal_plac) + sqrt(1-reliability_arousal)*rnorm(num) 
arousal_caff <- sqrt(reliability_arousal) * (true_arousal_caff) + sqrt(1-reliability_arousal)*rnorm(num) 
performance_plac <- (1/(1+exp(-true_arousal_plac)))*(1/(1+exp(true_arousal_plac)))   #inverted u function of arousal
performance_caff <- (1/(1+exp(-true_arousal_caff)))*(1/(1+exp(true_arousal_caff)))   #inverted u function of arousal

						#experimental conditions are block randomized
reward <- rep(c(0,1),num/2)   #reward vector of 0 or 1
punish <- rep(c(0,0,1,1),num/4) #punishment vector
block <- sort(rep(seq(1:(num/4)),4))  #block the data to allow for block randomization
temp.condition <- data.frame(reward,punish,block,rand =rnorm(num))  #experimental conditions ordered by block
condition <- temp.condition[order(temp.condition$block,temp.condition$rand),1:2]  #conditions are now block randomized


#true affect measures are a function of a trait, a situation, and their interaction 
TruePosAffect <- 1/(1+exp(-(weight_e * true_e + weightreward * condition$reward +weight_er * true_e * condition$reward)))
TrueNegAffect <- 1/(1+exp(-(weight_n * true_n + weightpunish * condition$punish +weight_np * true_n * condition$punish )))
TruePosAffect <- scale(TruePosAffect)     #standardize TruePosAffect to put on the same metric as the error scores 
TrueNegAffect <- scale(TrueNegAffect)     #standardize TrueNegAffect to put on the same metric as the error scores

#observed affect is a function of true affect and errors in measurement
PosAffect <- sqrt(reliability_P) * TruePosAffect+sqrt(1-reliability_P)*rnorm(num) 
NegAffect <- sqrt(reliability_N) * TrueNegAffect+sqrt(1-reliability_N)*rnorm(num) 

Data Organization

The data have been simulated and are ready for analysis. When doing the following analyses it is simpler if the data are organized in "data.frames" which are analogous to matrices except that different types of variables can be in different columns.

The analyses will be done as moderated multiple regression. This requires certain modifications to the data before putting them into a data frame. The most important is to center the Person Variable predictors around their mean and to convert the Experimental Variables to "factors" or categorical variables. The centering operation makes the product terms independent of the main effects.

In addition, some of the data should be transformed to allow for alernative demonstrations of analytic techniques. One transform is merely finding the sums and differences of the affect measures. This will be useful when we show that mixed within-between subject ANOVA can be thought of as a combination of two separate ANOVAS, one between subjects and one within subjects.



#organize all the data in a data frame to allow for analysis
#because it is also possible to do repeated measures as ANOVA on sums and differences
#the between effects are found by the sums
#the within effects are found the differences

PosplusNeg = PosAffect+NegAffect
PosminusNeg <- PosAffect - NegAffect
affect.data <- data.frame(extraversion,neuroticism,PosAffect,NegAffect,PosplusNeg,PosminusNeg,true_e,true_n,TruePosAffect,TrueNegAffect,
reward = as.factor(condition$reward),punish=as.factor(condition$punish))
centered.affect.data <- data.frame(scale(affect.data[,1:10],scale=FALSE),reward = as.factor(condition$reward),punish=as.factor(condition$punish))
drug.data <- data.frame(extraversion,true_arousal_plac,true_arousal_caff,arousal_plac,arousal_caff,performance_plac,performance_caff)


rm(reward,punish,block,condition,TruePosAffect,TrueNegAffect,PosAffect,NegAffect,extraversion, neuroticism, true_e, true_n,PosplusNeg,PosminusNeg,
true_arousal_plac,true_arousal_caff,arousal_plac,arousal_caff,performance_plac,performance_caff)
    #these have all been put into  data.frames and we don't need them individually anymore

The artificial data sets have been simulated and stored into three "data.frames" for easier manipulation. We can show the first 5 lines of the data.frame:

 print(affect.data[1:5,],digits=2)
  extraversion neuroticism PosAffect NegAffect PosplusNeg PosminusNeg true_e true_n TruePosAffect TrueNegAffect reward punish
1          4.9         4.5      1.63      2.45       4.08       -0.83    4.4    4.2          1.33           1.6      1      1
2          2.1         4.3     -0.95      1.56       0.61       -2.50    2.4    4.0         -0.97           1.5      0      1
3          3.3         2.1     -0.76     -0.91      -1.67        0.15    3.4    2.0         -0.97          -1.5      0      0
4          3.6         4.7      0.74      1.16       1.90       -0.43    3.6    4.8          1.17           1.1      1      0
5          3.0         2.1      0.38     -1.47      -1.09        1.85    3.4    2.3         -0.97          -1.1      0      0

The centered data have had the means removed:

 
 print(centered.affect.data[1:5,],digits=2)
  extraversion neuroticism PosAffect NegAffect PosplusNeg PosminusNeg true_e true_n TruePosAffect TrueNegAffect reward punish
1        1.916        1.61      1.58      2.51       4.09       -0.93   1.34   1.29          1.33           1.6      1      1
2       -0.911        1.42     -0.99      1.61       0.62       -2.61  -0.60   1.13         -0.97           1.5      0      1
3        0.372       -0.82     -0.80     -0.85      -1.66        0.05   0.33  -0.92         -0.97          -1.5      0      0
4        0.594        1.79      0.69      1.22       1.91       -0.53   0.60   1.94          1.17           1.1      1      0
5        0.059       -0.84      0.34     -1.41      -1.08        1.75   0.37  -0.58         -0.97          -1.1      0      0
> 

Basic descriptive statistics (using describe)

Basic descriptive statistics can be found by using the summary (), describe() and pairs.panels functions. (Describe and pairs.panels are part of the "psych" package.)

describe(affect.data)
              var   n  mean   sd median   min  max range   se
extraversion    1 100  2.96 1.07   3.05  0.33 5.50  5.16 0.11
neuroticism     2 100  2.92 0.85   2.88  0.92 5.26  4.34 0.09
PosAffect       3 100  0.05 1.09  -0.03 -2.43 2.70  5.13 0.11
NegAffect       4 100 -0.06 0.98   0.05 -2.11 2.45  4.56 0.10
PosplusNeg      5 100 -0.01 1.48   0.12 -2.73 4.08  6.81 0.15
PosminusNeg     6 100  0.10 1.46   0.04 -3.24 3.13  6.37 0.15
true_e          7 100  3.03 1.04   3.09  0.01 5.29  5.28 0.10
true_n          8 100  2.91 0.90   2.93  0.98 5.70  4.73 0.09
TruePosAffect   9 100  0.00 1.00  -0.48 -0.97 1.48  2.45 0.10
TrueNegAffect  10 100  0.00 1.00   0.26 -2.66 1.95  4.60 0.10
reward         11  NA    NA 0.50     NA    NA   NA    NA   NA
punish         12  NA    NA 0.50     NA    NA   NA    NA   NA


describe(drug.data)


                  var   n  mean   sd median   min  max range   se
extraversion        1 100  2.96 1.07   3.05  0.33 5.50  5.16 0.11
true_arousal_plac   2 100 -0.53 1.33  -0.45 -3.84 2.84  6.68 0.13
true_arousal_caff   3 100  0.51 1.10   0.51 -2.19 3.06  5.25 0.11
arousal_plac        4 100 -0.52 1.28  -0.43 -4.11 2.96  7.08 0.13
arousal_caff        5 100  0.38 1.04   0.36 -2.09 3.25  5.33 0.10
performance_plac    6 100  0.18 0.07   0.21  0.02 0.25  0.23 0.01
performance_caff    7 100  0.19 0.06   0.21  0.04 0.25  0.21 0.01

Basic descriptive statistics -- graphical presentation (using pairs.panels)

The data can also be visualized using the pairs.panels function (adapted from the pairs function of base R). splot plot of drug data

Part 3: Analyse the data using ANOVA and the general linear model.

Some of the analyses are shown merely to allow the user to appreciate the differences between various ways of model fitting.



#the first models do not use 0 centered data and are incorrect
#these are included merely to show what happens if the correct model is not used
#do the analyses using  a linear model without centering the data   ---- wrong
mod1w <- lm(PosAffect ~ extraversion+reward,data= affect.data)   #don't exam interactions
mod2w <- lm(NegAffect ~ neuroticism+punish,data = affect.data)
mod3w <- lm(PosAffect ~ extraversion*reward,data = affect.data) #look for interactions
mod4w <- lm(NegAffect ~ neuroticism*punish,data = affect.data)
mod5w <- lm(PosAffect ~ extraversion*neuroticism*reward*punish,data = affect.data)
				#show the results of the analyses
summary(mod1w,digits=2)
summary(mod2w,digits=2)
summary(mod3w,digits=2)
summary(mod4w,digits=2)
summary(mod5w,digits=2)


#do the analyses with centered data   -- this is the correct way -- note the differences 
#just look at main effects
mod1 <- lm(NegAffect ~ neuroticism+punish,data = centered.affect.data)   #just main effects
mod2 <- lm(PosAffect ~ extraversion+reward,data= centered.affect.data)   #don't exam interactions

#include interactions
# note that mod3 and mod4 are two different ways of specifying the interaction
mod3 <- lm(PosAffect ~ extraversion*reward,data = centered.affect.data) #look for interactions
mod4 <- lm(NegAffect ~ neuroticism+ punish + neuroticism*punish,data = centered.affect.data)
#go for the full models
mod5 <- lm(PosAffect ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)
mod6 <- lm(NegAffect ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)
mod6.5  <- lm(c(NegAffect,PosAffect) ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)


				#show the results of the analyses
summary(mod1,digits=2)
summary(mod2,digits=2)
summary(mod3,digits=2)
summary(mod4,digits=2)
summary(mod5,digits=2)
summary(mod6,digits=2)


Produces the following output

The first two analyses are just looking for main effects and are not as interesting as the third one.

 summary(mod1,digits=2)   #just the main effects of neuroticism and the negative film

Call:
lm(formula = NegAffect ~ neuroticism + punish, data = centered.affect.data)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.570842 -0.313535  0.004136  0.329371  1.269003 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.54682    0.07814  -6.998 3.37e-10 ***
neuroticism  0.72367    0.06514  11.110  < 2e-16 ***
punish1      1.09364    0.11053   9.895 2.26e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.5524 on 97 degrees of freedom
Multiple R-Squared: 0.6891,	Adjusted R-squared: 0.6827 
F-statistic: 107.5 on 2 and 97 DF,  p-value: < 2.2e-16 

> summary(mod2,digits=2)   #just extraversion + rewarding film

Call:
lm(formula = PosAffect ~ extraversion + reward, data = centered.affect.data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.13532 -0.46240  0.09547  0.43269  1.93528 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.84505    0.09716  -8.697 8.67e-14 ***
extraversion  0.13151    0.06432   2.045   0.0436 *  
reward1       1.69010    0.13745  12.296  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6869 on 97 degrees of freedom
Multiple R-Squared: 0.6133,	Adjusted R-squared: 0.6054 
F-statistic: 76.93 on 2 and 97 DF,  p-value: < 2.2e-16 

> summary(mod3,digits=2)   #this one includes the interaction, compare it to model 2

Call:
lm(formula = PosAffect ~ extraversion * reward, data = centered.affect.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0622 -0.4644  0.0832  0.4447  2.0437 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          -0.840116   0.095747  -8.774 6.37e-14 ***
extraversion         -0.005329   0.093516  -0.057   0.9547    
reward1               1.689352   0.135401  12.477  < 2e-16 ***
extraversion:reward1  0.252946   0.127146   1.989   0.0495 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6766 on 96 degrees of freedom
Multiple R-Squared: 0.6286,	Adjusted R-squared: 0.617 
F-statistic: 54.17 on 3 and 96 DF,  p-value: < 2.2e-16 

> summary(mod4,digits=2)   # the interaction is not significant -- compare to model 1

Call:
lm(formula = NegAffect ~ neuroticism + punish + neuroticism * 
    punish, data = centered.affect.data)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.565723 -0.309362  0.003068  0.330083  1.264219 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -0.54665    0.07855  -6.959 4.21e-10 ***
neuroticism          0.71679    0.08752   8.190 1.12e-12 ***
punish1              1.09369    0.11110   9.845 3.21e-16 ***
neuroticism:punish1  0.01561    0.13188   0.118    0.906    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.5552 on 96 degrees of freedom
Multiple R-Squared: 0.6892,	Adjusted R-squared: 0.6795 
F-statistic: 70.95 on 3 and 96 DF,  p-value: < 2.2e-16 


We can see this relationship more clearly if we graph it. In this plot, we show both the "true" data as well as the observed data. Plotting commands for this are given in the graphics section.

Figure 1: The effect of negative feedback or a sad movie upon Negative Affect

Figure 2: The effect of positive feedback or a happy movie upon Positive Affect

> summary(mod5,digits=2)

Call:
lm(formula = PosAffect ~ extraversion * neuroticism * reward * 
    punish, data = centered.affect.data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.12045 -0.41413  0.07548  0.48669  2.00568 

Coefficients:
                                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              -0.63697    0.14799  -4.304 4.51e-05 ***
extraversion                              0.07457    0.14006   0.532   0.5958    
neuroticism                              -0.04179    0.17649  -0.237   0.8134    
reward1                                   1.45328    0.20803   6.986 6.20e-10 ***
punish1                                  -0.42974    0.20259  -2.121   0.0369 *  
extraversion:neuroticism                 -0.04462    0.12356  -0.361   0.7189    
extraversion:reward1                      0.09537    0.20149   0.473   0.6372    
neuroticism:reward1                      -0.03706    0.23454  -0.158   0.8748    
extraversion:punish1                     -0.24191    0.21731  -1.113   0.2688    
neuroticism:punish1                       0.05500    0.28466   0.193   0.8473    
reward1:punish1                           0.44668    0.29069   1.537   0.1281    
extraversion:neuroticism:reward1         -0.19181    0.21128  -0.908   0.3666    
extraversion:neuroticism:punish1         -0.03263    0.29695  -0.110   0.9128    
extraversion:reward1:punish1              0.39800    0.28834   1.380   0.1711    
neuroticism:reward1:punish1              -0.12139    0.36464  -0.333   0.7400    
extraversion:neuroticism:reward1:punish1  0.27421    0.36118   0.759   0.4498    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.6821 on 84 degrees of freedom
Multiple R-Squared: 0.6698,	Adjusted R-squared: 0.6108 
F-statistic: 11.36 on 15 and 84 DF,  p-value: 1.618e-14 

> summary(mod6,digits=2)

Call:
lm(formula = NegAffect ~ extraversion * neuroticism * reward * 
    punish, data = centered.affect.data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.14498 -0.33986  0.01124  0.34204  1.20472 

Coefficients:
                                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)                              -0.75119    0.11463  -6.553 4.29e-09 ***
extraversion                              0.24492    0.10848   2.258  0.02656 *  
neuroticism                               0.79015    0.13670   5.780 1.23e-07 ***
reward1                                   0.44853    0.16113   2.784  0.00664 ** 
punish1                                   1.33237    0.15692   8.491 6.28e-13 ***
extraversion:neuroticism                  0.02834    0.09570   0.296  0.76787    
extraversion:reward1                     -0.16061    0.15606  -1.029  0.30638    
neuroticism:reward1                      -0.04336    0.18167  -0.239  0.81193    
extraversion:punish1                     -0.39715    0.16832  -2.360  0.02062 *  
neuroticism:punish1                      -0.04101    0.22048  -0.186  0.85290    
reward1:punish1                          -0.48703    0.22516  -2.163  0.03338 *  
extraversion:neuroticism:reward1         -0.27373    0.16365  -1.673  0.09811 .  
extraversion:neuroticism:punish1          0.06403    0.23000   0.278  0.78138    
extraversion:reward1:punish1              0.18717    0.22333   0.838  0.40437    
neuroticism:reward1:punish1               0.07575    0.28244   0.268  0.78919    
extraversion:neuroticism:reward1:punish1  0.20111    0.27975   0.719  0.47420    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.5283 on 84 degrees of freedom
Multiple R-Squared: 0.7538,	Adjusted R-squared: 0.7098 
F-statistic: 17.14 on 15 and 84 DF,  p-value: < 2.2e-16 



Note that the last two analyses could have been done in parallel by combining them into one model:


mod6.5  <- lm(c(NegAffect,PosAffect) ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)
summary(mod6.5) 

Multiple Dependent Measures- MANOVA

In the case of multiple dependent measures, there are a number of ways to treat the data. Multivariate Analysis of Variance (MANOVA) exams whether the various predictor variables have an overall effect upon all the dependent variables. If they do, then one can proceed to examine the separate dependent variables.

Traditionally, when examining several dependent variables (e.g., positive and negative affect), one might be interested in whether there is a difference in the effects depending upon which dependent variable is examined. That is, is there an interaction of dependent variable by some treatment?



#try MANOVA


mfit <- manova(cbind(PosAffect,NegAffect)~extraversion*neuroticism*reward*punish,data=centered.affect.data)
summary(mfit)
summary.aov(mfit)    #gives the "simple fits" for both dependent variables


The Manova output gives an overall test of the model, as well as separate estimates for each of the dependent variables:

mfit <- manova(cbind(PosAffect,NegAffect)~extraversion*neuroticism*reward*punish,data=centered.affect.data)
> summary(mfit)
                                       Df Pillai approx F num Df den Df    Pr(>F)    
extraversion                            1  0.037    1.612      2     83   0.20561    
neuroticism                             1  0.623   68.525      2     83 < 2.2e-16 ***
reward                                  1  0.644   75.235      2     83 < 2.2e-16 ***
punish                                  1  0.577   56.633      2     83 3.081e-16 ***
extraversion:neuroticism                1  0.012    0.525      2     83   0.59339    
extraversion:reward                     1  0.090    4.123      2     83   0.01963 *  
neuroticism:reward                      1  0.009    0.372      2     83   0.69070    
extraversion:punish                     1  0.084    3.801      2     83   0.02633 *  
neuroticism:punish                      1  0.001    0.022      2     83   0.97834    
reward:punish                           1  0.099    4.538      2     83   0.01348 *  
extraversion:neuroticism:reward         1  0.013    0.529      2     83   0.59089    
extraversion:neuroticism:punish         1  0.038    1.634      2     83   0.20141    
extraversion:reward:punish              1  0.023    0.961      2     83   0.38663    
neuroticism:reward:punish               1  0.004    0.152      2     83   0.85887    
extraversion:neuroticism:reward:punish  1  0.011    0.461      2     83   0.63195    
Residuals                              84                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> summary.aov(mfit)    #gives the "simple fits" for both dependent variables
 Response PosAffect :
                                       Df    Sum Sq   Mean Sq   F value  Pr(>F)    
extraversion                            1     1.254     1.254    2.6965 0.10431    
neuroticism                             1     0.718     0.718    1.5442 0.21746    
reward                                  1    70.663    70.663  151.8904 < 2e-16 ***
punish                                  1     0.659     0.659    1.4165 0.23733    
extraversion:neuroticism                1     0.342     0.342    0.7347 0.39381    
extraversion:reward                     1     2.531     2.531    5.4404 0.02207 *  
neuroticism:reward                      1     0.012     0.012    0.0263 0.87151    
extraversion:punish                     1     0.059     0.059    0.1270 0.72244    
neuroticism:punish                      1 5.367e-08 5.367e-08 1.154e-07 0.99973    
reward:punish                           1     1.321     1.321    2.8399 0.09566 .  
extraversion:neuroticism:reward         1     0.069     0.069    0.1486 0.70090    
extraversion:neuroticism:punish         1     0.465     0.465    0.9997 0.32024    
extraversion:reward:punish              1     0.779     0.779    1.6752 0.19912    
neuroticism:reward:punish               1     0.125     0.125    0.2687 0.60558    
extraversion:neuroticism:reward:punish  1     0.268     0.268    0.5764 0.44985    
Residuals                              84    39.079     0.465                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 Response NegAffect :
                                       Df Sum Sq Mean Sq  F value    Pr(>F)    
extraversion                            1  0.292   0.292   1.0449  0.309626    
neuroticism                             1 35.812  35.812 128.3103 < 2.2e-16 ***
reward                                  1  0.616   0.616   2.2080  0.141039    
punish                                  1 29.501  29.501 105.6980 < 2.2e-16 ***
extraversion:neuroticism                1  0.049   0.049   0.1750  0.676764    
extraversion:reward                     1  0.458   0.458   1.6413  0.203671    
neuroticism:reward                      1  0.184   0.184   0.6591  0.419156    
extraversion:punish                     1  2.144   2.144   7.6801  0.006874 ** 
neuroticism:punish                      1  0.012   0.012   0.0430  0.836194    
reward:punish                           1  1.344   1.344   4.8154  0.030968 *  
extraversion:neuroticism:reward         1  0.286   0.286   1.0251  0.314211    
extraversion:neuroticism:punish         1  0.776   0.776   2.7799  0.099181 .  
extraversion:reward:punish              1  0.150   0.150   0.5378  0.465375    
neuroticism:reward:punish               1  0.003   0.003   0.0117  0.914114    
extraversion:neuroticism:reward:punish  1  0.144   0.144   0.5168  0.474199    
Residuals                              84 23.445   0.279                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 



Multiple Dependent Measures- Repeated Measures

Traditionally, when examining several dependent variables (e.g., positive and negative affect), one might be interested in whether there is a difference in the effects depending upon which dependent variable is examined. That is, is there an interaction of dependent variable by some treatment? Until recently this has led to the use of repeated measures analysis of variance. The first two examples (models 7 and 8) demonstrate that the repeated measures can be broken up into two separate analyses, one looking at the sums of the scores, one looking at the differences of the scores. Compare these two to the repeated measures model (mod9). Although the Sums of Squares differ for the two models, note how the Mean Squares, F values, and resulting P values are identical for these two ways of thinking about the data.


#models 7 and 8 are included to show how there are alternative ways of doing a within subjects analysis 
mod7 <- lm(PosplusNeg  ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)    # this does the equivalent of the between part of a within-between analysis 
mod8 <- lm(PosminusNeg  ~ extraversion*neuroticism*reward*punish,data = centered.affect.data)    # this does the equivalent of the within part of a within-between analysis 

summary(mod7,digits=2)
summary(mod8,digits=2)

#repeated measures ANOVA 
#this requires a bit more work, in that we need to string out the data
attach(centered.affect.data)

PosNeg.df <- data.frame(PosAffect,NegAffect)   #make this into a 2 column data frame
stacked.pn <- stack(PosNeg.df)                 #string the data out for analysis
design.df <- data.frame(stacked.pn,rbind(centered.affect.data,centered.affect.data),subj=gl(num,1,2*num))
#note that subject is a factor in the repeated measures design
#the gl() function creates the factors 
#design.df        #show the data frame to understand it
detach(centered.affect.data)


mod7 <- aov(values ~ ind*extraversion*neuroticism*reward*punish+Error(subj),data=design.df)
summary(mod7)

#also possible to do repeated measures as ANOVA on sums and differences
#the between effects are found by the sums
#the within effects are found the differences
#thus, the following two analyses are the same as model 7
mod8 <- aov(PosplusNeg ~extraversion*neuroticism*reward*punish,data=centered.affect.data)
#
mod9 <- aov(PosminusNeg ~extraversion*neuroticism*reward*punish,data=centered.affect.data)
summary(mod8)    #the between part of model 7
summary(mod9)    #the within part of model 7


These analyses produce the following output:


summary(mod7)

Error: subj
                                       Df    Sum Sq   Mean Sq F value    Pr(>F)    
extraversion                            1     1.378     1.378  3.1773   0.07828 .  
neuroticism                             1    13.193    13.193 30.4219 3.763e-07 ***
reward                                  1    42.239    42.239 97.3986 1.063e-15 ***
punish                                  1    10.671    10.671 24.6055 3.623e-06 ***
extraversion:neuroticism                1     0.066     0.066  0.1524   0.69720    
extraversion:reward                     1     0.418     0.418  0.9633   0.32916    
neuroticism:reward                      1     0.051     0.051  0.1168   0.73341    
extraversion:punish                     1     1.457     1.457  3.3602   0.07033 .  
neuroticism:punish                      1     0.006     0.006  0.0138   0.90682    
reward:punish                           1 4.866e-05 4.866e-05  0.0001   0.99157    
extraversion:neuroticism:reward         1     0.318     0.318  0.7338   0.39408    
extraversion:neuroticism:punish         1     1.221     1.221  2.8160   0.09705 .  
extraversion:reward:punish              1     0.807     0.807  1.8603   0.17624    
neuroticism:reward:punish               1     0.044     0.044  0.1013   0.75108    
extraversion:neuroticism:reward:punish  1     0.403     0.403  0.9290   0.33789    
Residuals                              84    36.428     0.434                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Within
                                           Df    Sum Sq   Mean Sq   F value    Pr(>F)    
ind                                         1 3.755e-31 3.755e-31 1.209e-30  1.000000    
ind:extraversion                            1    0.1682    0.1682    0.5415  0.463880    
ind:neuroticism                             1   23.3372   23.3372   75.1217 2.774e-13 ***
ind:reward                                  1   29.0406   29.0406   93.4806 2.683e-15 ***
ind:punish                                  1   19.4891   19.4891   62.7346 8.747e-12 ***
ind:extraversion:neuroticism                1    0.3245    0.3245    1.0446  0.309678    
ind:extraversion:reward                     1    2.5713    2.5713    8.2770  0.005089 ** 
ind:neuroticism:reward                      1    0.1456    0.1456    0.4686  0.495521    
ind:extraversion:punish                     1    0.7454    0.7454    2.3995  0.125136    
ind:neuroticism:punish                      1    0.0060    0.0060    0.0194  0.889542    
ind:reward:punish                           1    2.6651    2.6651    8.5790  0.004377 ** 
ind:extraversion:neuroticism:reward         1    0.0370    0.0370    0.1191  0.730890    
ind:extraversion:neuroticism:punish         1    0.0198    0.0198    0.0636  0.801448    
ind:extraversion:reward:punish              1    0.1227    0.1227    0.3949  0.531428    
ind:neuroticism:reward:punish               1    0.0843    0.0843    0.2715  0.603714    
ind:extraversion:neuroticism:reward:punish  1    0.0095    0.0095    0.0307  0.861399    
Residuals                                  84   26.0954    0.3107                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> #also possible to do repeated measures as ANOVA on sums and differences
> #the between effects are found by the sums
> #the within effects are found the differences
> #thus, the following two analyses are the same as model 7
> mod8 <- aov(PosplusNeg ~extraversion*neuroticism*reward*punish,data=centered.affect.data)
> #
> mod9 <- aov(PosminusNeg ~extraversion*neuroticism*reward*punish,data=centered.affect.data)
> summary(mod8)    #the between part of model 7
                                       Df    Sum Sq   Mean Sq F value    Pr(>F)    
extraversion                            1     2.756     2.756  3.1773   0.07828 .  
neuroticism                             1    26.386    26.386 30.4219 3.763e-07 ***
reward                                  1    84.477    84.477 97.3986 1.063e-15 ***
punish                                  1    21.341    21.341 24.6055 3.623e-06 ***
extraversion:neuroticism                1     0.132     0.132  0.1524   0.69720    
extraversion:reward                     1     0.836     0.836  0.9633   0.32916    
neuroticism:reward                      1     0.101     0.101  0.1168   0.73341    
extraversion:punish                     1     2.914     2.914  3.3602   0.07033 .  
neuroticism:punish                      1     0.012     0.012  0.0138   0.90682    
reward:punish                           1 9.732e-05 9.732e-05  0.0001   0.99157    
extraversion:neuroticism:reward         1     0.636     0.636  0.7338   0.39408    
extraversion:neuroticism:punish         1     2.442     2.442  2.8160   0.09705 .  
extraversion:reward:punish              1     1.613     1.613  1.8603   0.17624    
neuroticism:reward:punish               1     0.088     0.088  0.1013   0.75108    
extraversion:neuroticism:reward:punish  1     0.806     0.806  0.9290   0.33789    
Residuals                              84    72.856     0.867                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> summary(mod9)    #the within part of model 7
                                       Df Sum Sq Mean Sq F value    Pr(>F)    
extraversion                            1  0.336   0.336  0.5415  0.463880    
neuroticism                             1 46.674  46.674 75.1217 2.774e-13 ***
reward                                  1 58.081  58.081 93.4806 2.683e-15 ***
punish                                  1 38.978  38.978 62.7346 8.747e-12 ***
extraversion:neuroticism                1  0.649   0.649  1.0446  0.309678    
extraversion:reward                     1  5.143   5.143  8.2770  0.005089 ** 
neuroticism:reward                      1  0.291   0.291  0.4686  0.495521    
extraversion:punish                     1  1.491   1.491  2.3995  0.125136    
neuroticism:punish                      1  0.012   0.012  0.0194  0.889542    
reward:punish                           1  5.330   5.330  8.5790  0.004377 ** 
extraversion:neuroticism:reward         1  0.074   0.074  0.1191  0.730890    
extraversion:neuroticism:punish         1  0.040   0.040  0.0636  0.801448    
extraversion:reward:punish              1  0.245   0.245  0.3949  0.531428    
neuroticism:reward:punish               1  0.169   0.169  0.2715  0.603714    
extraversion:neuroticism:reward:punish  1  0.019   0.019  0.0307  0.861399    
Residuals                              84 52.191   0.621                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 



The best way to show this set of interactions is by a four panel graph. (Details on making this graph are given in the graphics section.) Note the clear additive effects of Neuroticism and Negative movies on Negative Affect (top left panel), the lack of an effect of Extraversion on Negative Affect (top right), the effect of positive movies on Positive affect (bottom panels), and the effect of extraversion by movie condition (bottom right). There is an interaction between movie type (reward x punishment) and affect measure, shown in the within subjects analysis.

extraversion x neuroticism x movie condition

Study 4: The effect of caffeine and extraversion on cognitive performance

Analyse the drug by extraversion data. Since this is a within subjects design, we will need to rearrange the data slightly.


#repeated measures ANOVA 
#this requires a bit more work, in that we need to string out the data

attach(drug.data)

drug.df <- data.frame(performance_plac,performance_caff)   #make this into a 2 column data frame
stacked.drug <- stack(drug.df)                 #string the data out for analysis
drug.design.df <- data.frame(stacked.drug,rbind(drug.data,drug.data),subj=gl(num,1,2*num))
#note that subject is a factor in the repeated measures design
#the gl() function creates the factors 
#drug.design.df        #show the data frame to understand it
detach(drug.data)

mod10 <- aov(values ~ ind*extraversion+Error(subj),data=drug.design.df)
summary(mod10)


summary(mod10)

Error: subj
             Df  Sum Sq Mean Sq F value Pr(>F)
extraversion  1 0.00506 0.00506  1.4492 0.2316
Residuals    98 0.34197 0.00349               

Error: Within
                 Df  Sum Sq Mean Sq F value   Pr(>F)   
ind               1 0.00568 0.00568   1.506 0.222685   
ind:extraversion  1 0.04068 0.04068  10.791 0.001415 **
Residuals        98 0.36940 0.00377                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

This can be shown graphically in two different ways. The left hand panel simply shows the two regression lines (Performance as a function of extraversion for Placebo and Caffeine conditions). This shows that performance declines with extraversion in the placebo condition and increases in the caffeine condition. This is the within subjects interaction with extraversion term in the regression analysis. The right hand panel clarifies this relationship by using a "lowess" smooth of the data, suggesting a curvilinear fit for both drug conditions with a difference of the peaks of the curves.

extraversion by drug condition







Graphical Data Summaries

Graphical presentations of the data make these results easier to understand. With not many subjects, it is useful to show the data points as well, but as the number of participants increases, it is probably better to show just the regression lines. Here we show the commands that made the graphs presented above.

For an overview of how to use R for graphics, consult the short guide to R, for detailed discussions of graphics in R, consult the many help pages for graphics.



#graphical output helps understand interactions
#the following code shows how to draw these graphs
#First, Figure 1

                               
symb <- c(19,25,3,23)                              #choose some nice plotting symbols
colors <- c("black","red","green","blue")          #choose some nice colors

par(mfrow=c(1,2))     #specify the number of rows and columns in a multi panel graph  



#make a two panel graph of the neuroticism x punishment study

attach(affect.data) 	#this gives us direct access to the variables
#first draw the latent results, fit them, and label them 
plot(true_n,TrueNegAffect,pch=symb[punish],col=colors[punish],bg=colors[punish],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Negative Affect",
xlab="Latent Neuroticism",main="True relationships")               #show the data points
by(affect.data,punish,function(x) abline(lm(TrueNegAffect~true_n,data=x)))
text(3,-1.5,"Neutral")
text(3,2,"Negative")
by(affect.data,punish,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2))

plot(neuroticism,NegAffect,pch=symb[punish],col=colors[punish],bg=colors[punish],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Negative Affect",
xlab="Observed Neuroticism", main="Observed relationships")               #show the data points
by(affect.data,punish,function(x) abline(lm(NegAffect~neuroticism,data=x)))
text(3,-2.3,"Neutral")
text(3,2,"Negative")
by(affect.data,punish,function(x) summary(lm(NegAffect~neuroticism,data=x),digits=2))
detach(affect.data)                 #otherwise, we have copies floating around



Figure 1: The effect of negative feedback or a sad movie upon Negative Affect #now do a 2 panel graph of the extraversion x reward study #first draw the latent results, then the "real results" attach(affect.data) #this gives us direct access to the variables colors <- c("black","blue","green","red") plot(true_e,TruePosAffect,pch=symb[reward],col=colors[reward],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Positive Affect", xlab="Latent extraversion",main="True relationships") #show the data points by(affect.data,reward,function(x) abline(lm(TruePosAffect~true_e,data=x))) text(3,-1.5,"Neutral") text(3,2,"Positive") by(affect.data,reward,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2)) plot(extraversion,PosAffect,pch=symb[reward],col=colors[reward],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Positive Affect", xlab="Observed Extraversion", main="Observed relationships") #show the data points by(affect.data,reward,function(x) abline(lm(PosAffect~extraversion,data=x))) text(3,-2.3,"Neutral") text(3,2,"Positive") by(affect.data,reward,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2)) detach(affect.data) #otherwise, we have copies floating around

Figure 2: The effect of positive feedback or a happy movie upon Positive Affect

Now, draw a 2 x 2 panel graph of the complete relationships. First we will do this for the latent (pure) data, and then for the observed data.



#graphical output helps understand interactions

attach(affect.data) 	#this gives us direct access to the variables
                               
symb <- c(19,25,3,23)                              #choose some nice plotting symbols
colors <- c("black","red","blue","green")          #choose some nice colors

par(mfrow=c(2,2))    
#make a four panel graph of the personality x film data 
#the rows will be affect
#the columns personality variable
#The four separate lines in each graph represent the four movie conditions
#solid lines compare neutral vs. negative (threat) movies
#dashed lines compare neutral vs. positive (reward) movies

moviecond <- as.numeric(punish) + 2* as.numeric(reward) -2   #movie conditions
#movie conditions 1 = neutral   2 = punish alone  3= reward alone   4 = punish+reward

#first draw the latent results, fit them, and label them 
plot(true_n,TrueNegAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[moviecond],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Negative Affect",
xlab="Latent Neuroticism",main="Neuroticism x Movies: Negative Affect")               #show  the data points

by(affect.data,moviecond,function(x) abline(lm(TrueNegAffect~true_n,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Negative")

plot(true_e,TrueNegAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[punish],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Negative Affect",
xlab="Latent Extraversion",main="Extraversion x Movies: Negative Affect")               #show  the data points
by(affect.data,moviecond,function(x) abline(lm(TrueNegAffect~true_e,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Negative")

plot(true_n,TruePosAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Positive Affect",
xlab="Latent Neuroticism",main="Neuroticism x Movies: Positive Affect")               #show  the data points
by(affect.data,moviecond,function(x) abline(lm(TruePosAffect~true_n,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Positive")

plot(true_e,TruePosAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Positive Affect",
xlab="Latent Extraversion",main="Extraversion x Movies: Positive Affect")               #show  the data points
by(affect.data,moviecond,function(x) abline(lm(TruePosAffect~true_e,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Positive")

The latent data (without error) show the effects of movie conditions upon positive and negative affect.

extraversion x neuroticism x pa x na latent

#now do it for the real data

symb <- c(19,25,3,23)                              #choose some nice plotting symbols
colors <- c("black","red","blue","green")          #choose some nice colors
attach(affect.data)

par(mfrow=c(2,2))    
#make a four panel graph of the personality x film data 
#the rows will be affect
#the columns personality variable
#The four separate lines in each graph represent the four movie conditions
#solid lines compare neutral vs. negative (threat) movies
#dashed lines compare neutral vs. positive (reward) movies

moviecond <- as.numeric(punish) + 2* as.numeric(reward) -2   #movie conditions
#movie conditions 1 = neutral   2 = punish alone  3= reward alone   4 = punish+reward

#first draw the latent results, fit them, and label them 
plot(neuroticism,NegAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[moviecond],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Negative Affect",
xlab="Observed Neuroticism",main="Neuroticism x Movies: Negative Affect")               #show  the data points

by(affect.data,punish,function(x)      abline(lm(NegAffect~neuroticism,data=x),lty="dashed"))  #add regression slopes

text(3,-1.5,"Neutral")
text(3,2,"Negative")

plot(extraversion,NegAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[punish],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Negative Affect",
xlab="Observed Extraversion",main="Extraversion x Movies: Negative Affect")               #show  the data points
by(affect.data,punish,function(x) abline(lm(NegAffect~extraversion,data=x),lty="dashed"))

text(3,-1.5,"Neutral")
text(3,2,"Negative")


plot(neuroticism,PosAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Positive Affect",
xlab="Observed Neuroticism",main="Neuroticism x Movies: Positive Affect")               #show  the data points
by(affect.data,reward,function(x) abline(lm(PosAffect~neuroticism,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Positive")

plot(extraversion,PosAffect,pch=symb[moviecond],col=colors[moviecond],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Positive Affect",
xlab="Observed Extraversion",main="Extraversion x Movies: Positive Affect")               #show  the data points
by(affect.data,reward,function(x) abline(lm(PosAffect~extraversion,data=x)))

text(3,-1.5,"Neutral")
text(3,2,"Positive")

detach(affect.data)                 #otherwise, we have copies floating around



by(affect.data,punish,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2))

plot(neuroticism,NegAffect,pch=symb[punish],col=colors[punish],bg=colors[punish],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Negative Affect",
xlab="Observed Neuroticism", main="Observed relationships")               #show  the data points
by(affect.data,punish,function(x) abline(lm(NegAffect~neuroticism,data=x)))
text(3,-2.3,"Neutral")
text(3,2,"Negative")
by(affect.data,punish,function(x) 
summary(lm(NegAffect~neuroticism,data=x),digits=2))
detach(affect.data)                 #otherwise, we have copies floating around



#now do a 2 panel graph of the extraversion x reward study
#first draw the latent results, then the "real results"

attach(affect.data) 	#this gives us direct access to the variables
colors <- c("black","blue","green","red")  
plot(true_e,TruePosAffect,pch=symb[reward],col=colors[reward],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Latent Positive Affect",
xlab="Latent extraversion",main="True relationships")               #show  the data points
by(affect.data,reward,function(x) abline(lm(TruePosAffect~true_e,data=x)))
text(3,-1.5,"Neutral")
text(3,2,"Positive")
by(affect.data,reward,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2))

plot(extraversion,PosAffect,pch=symb[reward],col=colors[reward],bg=colors[reward],cex=1.0,ylim =c(-2.5,2.5),ylab="Observed Positive Affect",
xlab="Observed Extraversion", main="Observed relationships")               #show  the data points
by(affect.data,reward,function(x) abline(lm(PosAffect~extraversion,data=x)))
text(3,-2.3,"Neutral")
text(3,2,"Positive")
by(affect.data,reward,function(x) summary(lm(PosAffect~extraversion,data=x),digits=2))
detach(affect.data)                 #otherwise, we have copies floating around


The "real" data are much noiser than the latent data but the same pattern of interactions is still shown.







Show the extraversion x caffeine graph. Note how we are fitting both linear and lowess fits to the data. Which provides a better description?



#graph the caffeine and placebo by extraversion data

attach(drug.design.df)


par(mfrow=c(1,2))                            #set the number of rows and columns
symb <- c(19,25,3,23)                              #choose some nice plotting symbols
colors <- c("black","red","blue","green")          #choose some nice colors

plot(extraversion,values,pch=symb[ind],col=colors[ind],cex=1.0,ylim =c(-0,.3),ylab="Performance",xlab=" Extraversion",main="Linear fits")      
by(drug.design.df,ind,function(x) abline(lm(values~extraversion,data=x)))
text(5,.16,"Placebo")
text(5,.20,"Caffeine")

plot(extraversion,values,pch=symb[ind],col=colors[ind],cex=1.0,ylim =c(-0,.3),ylab="Performance",xlab=" Extraversion",main="Lowess fits")      

text(5,.16,"Placebo")
text(5,.2,"Caffeine")
lines(lowess(performance_plac~extraversion))
lines(lowess(performance_caff~extraversion))

detach(drug.design.df)   #otherwise we have copies floating around
extraversion by drug condition







Multiple Dependent Measures- Multilevel or Mixed Model approaches

Traditionally, when examining several dependent variables (e.g., positive and negative affect), one might be interested in whether there is a difference in the effects depending upon which dependent variable is examined. That is, is there an interaction of dependent variable by some treatment? Although traditionally repeated measures ANOVA has been used (see above), more recent treatments are using a variety of procedures that generically are studying "mixed effects". Known as HLM, Proc Mixed, or just nmle, these procedures attempt to model the data at multiple levels, the subject level (where the repeated measures occur) and at a higher level (other characteristics of the subject.) See Pinheiro, J. C. & Bates, D.M. (2000) Mixed-effects models in S and S-Plus. Springer, New York. for more detail.


##Newer developments allow for mixed model estimates

library(nmle)       #we need to use the nmle package

mod.lme <- lme(values ~ind*reward*punish, random = ~1 |subj,data=design.df)
mod.lme <- lme(values ~ind*reward*punish, random = ~1 |neuroticism/subj,data=design.df)
mod.lme <- lme(values ~ind*reward*punish*extraversion*neuroticism , random = ~1 |subj,data=design.df)



Simulation 3: Challenges in measurement

A serious problem in inference about main effects and interactions is the relationship between observed variables and latent variables. If the OV is a monotonic but non linear function of the latent variable, and the latent variable is a function of some latent trait as well as some experimental manipulation, then it is likely that some observed interactions are due to the non linearity of the function rather than a true interaction.

In the chapter we discuss the problem of inference and scaling artifact. The R code for drawing the graph presented there is shown here.

#use curve to draw theoretical curves
#use margin control to make multi-panel graphs
#note that we scale the entire figure from 0 to 100 


#pdf("revelle.figure3.pdf",height=6,width=9)

par(mfrow=c(1,3))
par(mar=c(5,2,4,0))
curve(100/(1+exp(-x+2)),0,1,ylim=c(0,100),xaxp=c(0,1,1),xlab="Time during Year",main="Reading", cex=1.5)
curve(100/(1+exp(-x+3)),0,1,add=TRUE)
curve(100/(1+exp(-x+1)),0,1,add=TRUE)
text(.5,5,"Junior College")
text(.5,15,"Non selective")
text(.5,35,"Selective")


par(mar=c(5,0,4,0))
curve(100/(1+exp(-x-2)),-0,1,ylim=c(0,100),yaxt="n",xaxp=c(0,1,1),xlab="Time during Year",main="Mathematics")
curve(100/(1+exp(-x-3)),-0,1,add=TRUE)
curve(100/(1+exp(-x-1)),-0,1,add=TRUE)
text(.5,75,"Junior College")
text(.5,90,"Non selective")
text(.5,100,"Selective")



par(mar=c(5,0,4,1))
curve(100/(1+exp(-x-2)),-3,3,ylim=c(0,100),xlim = c(-2.5,2.5),yaxt="n",xlab="latent growth",main="Observed by Latent",lty="dotted")
#curve(100/(1+exp(-x-3)),-3,3,add=TRUE,lty="dashed")
#curve(100/(1+exp(-x+3)),-3,3,add=TRUE,lty="dashed")
curve(100/(1+exp(-x+2)),-3,3,add=TRUE,lty="dotted")
#curve(100/(1+exp(-x-1)),-3,3,add=TRUE,lty="dashed")
#curve(100/(1+exp(-x+1)),-3,3,add=TRUE,lty="dashed")
curve(100/(1+exp(-x-2)),1,2,add=TRUE,lwd=4)
curve(100/(1+exp(-x+2)),1,2,add=TRUE,lwd=4)
curve(100/(1+exp(-x+2)),-0,1,add=TRUE) 
curve(100/(1+exp(-x-2)),-0,1,add=TRUE)  
curve(100/(1+exp(-x+2)),-1,0,add=TRUE,lty="dashed",lwd=2)
curve(100/(1+exp(-x-2)),-1,0,add=TRUE,lty="dashed",lwd=2)

text(.2,5,"Junior College")
text(1,15,"Non selective")
text(2,35,"Selective")
text(.2,80,"Junior College")
text(1,90,"Non selective")
text(1.8,95,"Selective")
#dev.off()

In the following simulations this point is made somewhat clearer. Following the work in Item Response Theory, the OV is a logistic function of a latent trait. Situational manipulations are simulated by shifts along the latent trait. Latent interactions are created as product terms.



#measurement issues are a serious problem in inference
#The following code generates raw data for latent true scores
#and the creates observed scores as logistic transforms of latent scores
#Both additive and interactive effects are shown

#first some arbitrary constants
nsubjects <- 200
additive <- -2
interaction <- 2
error <- 1
constant <- 0   #additive constant to make all variable positive

theta <- rnorm(nsubjects) + error*rnorm(nsubjects) + constant   #basic trait has error 

addtheta <- theta - additive
interact <- theta*interaction
observed1 <- 1/(1+exp(-theta + constant))
observed2 <- 1/(1+exp(additive - theta + constant))
observed3 <- 1/(1+exp(interaction*(additive -theta + constant)))
mydata <- data.frame(theta,addtheta,interact,observed1,observed2,observed3)
mydata.log <- log(mydata)
describe(mydata)
pairs.panels(mydata)



op <- par(mfrow=c(4,2))
plot(theta,theta,xlab="latent score", ylab="latent score",main="simple additive effects at latent level",ylim=c(-3,3))
points(theta,addtheta)
text(1,.5,"control")
text(-2.5,.5,"experimental")

plot(theta,theta ,xlab="latent score", ylab="latent score",main="simple interaction effects at latent level",ylim=c(-3,3))
points(theta,interact)
text(1,.5,"control")
text(0,1,"experimental")

plot(theta,observed1,xlab="latent score", ylab="observed score",main="observed scores with simple additive effects at latent level",ylim=c(0,1))
points(theta,observed2)
text(1,.5,"control")
text(-1.5,.5,"experimental")

plot(theta,observed1,xlab="latent score", ylab="observed score",main="observed scores simple interaction effects at latent level",ylim=c(0,1))
points(theta,observed3)
text(1,.5,"control")
text(-1.5,.5,"experimental")

plot(observed1,observed1,xlab="observed score", ylab="observed score",main="observed scores when latent have additive effects",ylim=c(0,1))
points(observed1,observed2)
text(.6,.5,"control")
text(.2,.5,"experimental")

plot(observed1,observed1,xlab="observed score", ylab="observed score",main="observed scores when latent have interaction effects",ylim=c(0,1))
points(observed1,observed3)
text(.6,.5,"control")
text(.2,.5,"experimental")



plot(log(observed1),log(observed1),xlab="log observed score", ylab="log observed score",main="log observed scores when latent have additive effects")
points(log(observed1),log(observed2))
text(.6,.5,"control")
text(.2,.5,"experimental")

plot(log(observed1),log(observed1),xlab="log observed score", ylab="log observed score",main="log observed scores when latent have interaction effects")
points(log(observed1),log(observed3))
text(.6,.5,"control")
text(.2,.5,"experimental")










logistic <- function(x,a=0,b=1) {1/(1+exp(b*(a-x)))}
raw <- function(x) {x}

op <- par(mfrow=c(3,2))

curve(raw(x),-3,3,101,ylab ="latent score",xlab="latent score",lty="dashed")
curve(raw(x+1),-3,3,101,add=T)
#curve(raw(x-1),-3,3,101,add=T)
text(1,.5,"control")
text(-2,.5,"experimental")
title("Simple additive effects")

curve(raw(x),-3,3,101,ylab ="latent score",xlab="latent score",lty="dashed")
#curve(raw(x*2),-3,3,101,add=T)
curve(raw((x+1)*2),-3,3,101,add=T)
text(1,.5,"control")
text(-2,.5,"experimental")
title("Interaction + additive effects")


curve(logistic(x),-3,3,101,xlab="latent score", ylab="observed score",lty="dashed")
curve(logistic(x,a=-1),-3,3,101,add=T)
#curve(logistic(x,a=-1),-3,3,101,add=T)
text(1,.5,"control")
text(-2,.5,"experimental")
title("Simple additive effects on the latent variable")

curve(logistic(x), -3,3,101,xlab="latent score", ylab="observed score",lty="dashed")
#curve(logistic(x,b=2), -3,3,101,add=T)
curve(logistic(x,a=-1,b=2),-3,3,101,add=T)
text(1,.5,"control")
text(-2,.5,"experimental")
title("Interactons +  additive effects on the latent variable")

curve(logistic(x),-4,0,101,xlab="latent score", ylab="observed score",lty="dashed",ylim=c(0,1))
curve(logistic(x,a=-2),-4,0,101,add=T)
#curve(logistic(x,a=-1),-4,0,101,add=T)
text(-1,.4,"control")
text(-3,.4,"experimental")
title("Simple additive effects on the latent variable")

curve(logistic(x),0,4,101,xlab="latent score", ylab="observed score",lty="dashed",ylim=c(0,1))
curve(logistic(x,a=-2),0,4,101,add=T)
#curve(logistic(x,a=-1),0,4,101,add=T)
text(3,.8,"control")
text(1,.8,"experimental")
title("Simple additive effects on the latent variable")
mtext("The effect of scaling on inference",3,outer=TRUE,line=1,cex=1.5) 




The reader is encouraged to take the code in these simulations and to generate more cases, with different parameters to see what happens.
part of a guide to R
Version of July 10, 2006
William Revelle
Department of Psychology
Northwestern University