Simulation as a tool for personality research
William Revelle
Northwestern University
To appear in: Personality Research Methods, edited by Robins, Fraley and Krueger
Contents:
- Using the R computer package as a tool for personality research
- Data simulation as a research tool
- Data organization and simple summary statistics
- Using the general linear model to fit data
- Graphical presentations techniques
- 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).
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.
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.
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.
#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
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