
R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

#there is a handy library called schoRsch that gives you partial eta squared from ezANOVA
#http://finzi.psych.upenn.edu/library/schoRsch/html/anova_out.html


> library(car)
Loading required package: carData
Warning messages:
1: package ‘car’ was built under R version 3.6.3 
2: package ‘carData’ was built under R version 3.6.3 
> library(reshape)
Warning message:
package ‘reshape’ was built under R version 3.6.3 
> 
> zz <- file.path("U:","My Documents")
> 
> File.Location <- file.path(zz,"rep_meas_formattednogroup.csv")
> y <- read.csv(File.Location)
> y <- data.frame(y)
> 
> # obtain mean and sd for y1 in each group
> 
> library(data.table)
data.table 1.13.0 using 3 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: ‘data.table’

The following object is masked from ‘package:reshape’:

    melt

Warning message:
package ‘data.table’ was built under R version 3.6.3 
> dt <- data.table(y)
> dt[,list(mean=mean(y1),sd=sd(y1)),by=group]
Error in eval(bysub, parent.frame(), parent.frame()) : 
  object 'group' not found
> 
> subject <- c(1:18)
> 
> # create a data frame which contains only the response and subject indicator
> isr <- data.frame(y,subject)
> 
> isr <- data.frame(y,subject)
> 
> # group means for y1; note: could gave created group within R
> y$group = rep(1:2,each=9)
> tapply(y$y1,y$group,mean)
       1        2 
1.888889 2.333333 
> tapply(y$y1,y$group,sd)
       1        2 
0.781736 1.000000 

> 
> 
> # can look at distributions of repeated measures (normaliity assumptions about the distribution of differences)
> # see later for  a method for assessing residuals
> boxplot(y$y2-y$y1)
> boxplot(y$y2-y$y1,subset=y$group==1)
> 
> # we will just put the response into long format so removing group from the data frame; will generate this again after we have put the response
> # into long format
> 
> # we need to get rid of the between subjects factor which we will recreate later for the formatting to work
> 
> isr$group <- NULL
> 
> library(reshape)
> mdata = melt(isr, id=c("subject"))
Warning message:
In melt(isr, id = c("subject")) :
  The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(isr). In the next version, this warning will become an error.
> colnames(mdata) = c("subject", "var", "score")
> mdata <- data.frame(mdata)
> #mdata <- arrange(mdata,subject)
> head(mdata)
  subject var score
1       1  y1     3
2       2  y1     2
3       3  y1     1
4       4  y1     1
5       5  y1     2
6       6  y1     3
> 
> 
> # more generally could construct the long formatted group  variable from the group in the repmeas_formatted.csv input file
> group2 <- rep(y$group,6)
> 
> mood <- rep(c(1:3),each=18)
> mood <- rep(mood, 2)
> time <- rep(c(1:2),each=54)
> 
> 
> #simple effect comparing mood (Within subjects factor) in group1 (Between subjects factor)
> 
> # Paired t-tests for three mood groups using all data
> mdata <- data.frame(mdata,mood,time)
> attach(mdata)
The following objects are masked _by_ .GlobalEnv:

    mood, subject, time

> 
> #install.packages("ez")
> library(ez)
Registered S3 methods overwritten by 'lme4':
  method                          from
  cooks.distance.influence.merMod car 
  influence.merMod                car 
  dfbeta.influence.merMod         car 
  dfbetas.influence.merMod        car 
Warning message:
package ‘ez’ was built under R version 3.6.3 
> 
> # there are several ways to do mixed ANOVA in R; Field recommends ezANOVA which is in the ez library
> # this produces the Greenhouse-Geisser corrected p-values
> 
> # We have to make sure that ezANOVA will recognize the group inputs as factors which we can do using recode
> # otherwise we get incorrect dfs and also would not obtain the Greenhouse-Geisser output
> 
> # factor usually works at constructing factors but doesn't seem to work with ezANOVA
> 
> #mdata$mood2 <- factor(mdata$mood, labels = c('A','B','C'))
> #mdata$time <- factor(mdata$time, labels=c('base','follow-up'))
> #mdata$group2 <- factor(mdata$group, labels=c('controls','patients'))
> 
> #  so we use recode which does
> 
> mdata$mood2 <- recode(mdata$mood, "1='A'; 2='B'; 3='C'")
> mdata$time <- recode(mdata$time, "1='base';2='follow-up'")
> #mdata$group2 <- recode(mdata$group, "1='controls';2='patients'")
> mdata$group2 <- recode(mdata$group2, "1='controls';2='patients'")
Error in `$<-.data.frame`(`*tmp*`, group2, value = numeric(0)) : 
  replacement has 0 rows, data has 108
> 
> 
> # can fit WW model
> out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),detailed=T,type=3)
Warning: Converting "subject" to factor for ANOVA.
Warning: Converting "mood2" to factor for ANOVA.
Warning: Converting "time" to factor for ANOVA.
> out
$ANOVA
       Effect DFn DFd         SSn      SSd           F            p p<.05
1 (Intercept)   1  17 1020.592593 31.40741 552.4198113 2.107911e-14     *
2       mood2   2  34   19.907407 46.09259   7.3423061 2.236055e-03     *
3        time   1  17   23.148148 37.51852  10.4886476 4.828863e-03     *
4  mood2:time   2  34    2.907407 62.42593   0.7917532 4.612252e-01      
         ges
1 0.85188735
2 0.10087267
3 0.11539882
4 0.01612075

$`Mauchly's Test for Sphericity`
      Effect         W          p p<.05
2      mood2 0.9842713 0.88088393      
4 mood2:time 0.7039040 0.06027085      

$`Sphericity Corrections`
      Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
2      mood2 0.9845149 0.002369657         * 1.1124543 0.002236055         *
4 mood2:time 0.7715478 0.433440137           0.8337931 0.441793866          

> zz <- file.path("U:","My Documents")
> 
> File.Location <- file.path(zz,"rep_meas_formattednogroup.csv")
> y <- read.csv(File.Location)
> y <- data.frame(y)
> 
> # obtain mean and sd for y1 in each group
> 
> library(data.table)
> dt <- data.table(y)
> dt[,list(mean=mean(y1),sd=sd(y1)),by=group]
Error in eval(bysub, parent.frame(), parent.frame()) : 
  object 'group' not found
> 
> subject <- c(1:18)
> 
> # create a data frame which contains only the response and subject indicator
> isr <- data.frame(y,subject)
> 
> library(reshape)
> mdata = melt(isr, id=c("subject"))
Warning message:
In melt(isr, id = c("subject")) :
  The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(isr). In the next version, this warning will become an error.
> colnames(mdata) = c("subject", "var", "score")
> mdata <- data.frame(mdata)
> #mdata <- arrange(mdata,subject)
> head(mdata)
  subject var score
1       1  y1     3
2       2  y1     2
3       3  y1     1
4       4  y1     1
5       5  y1     2
6       6  y1     3
> 
> 
> # more generally could construct the long formatted group  variable from the group in the repmeas_formatted.csv input file
> group2 <- rep(c(1,2),each=9)
> 
> mood <- rep(c(1:3),each=18)
> mood <- rep(mood, 2)
> time <- rep(c(1:2),each=54)
> 
> 
> #simple effect comparing mood (Within subjects factor) in group1 (Between subjects factor)
> 
> # Paired t-tests for three mood groups using all data
> mdata <- data.frame(mdata,mood,time,group2)
> attach(mdata)
The following objects are masked _by_ .GlobalEnv:

    group2, mood, subject, time

The following objects are masked from mdata (pos = 4):

    mood, score, subject, time, var

> 
> #install.packages("ez")
> library(ez)
> 
> # there are several ways to do mixed ANOVA in R; Field recommends ezANOVA which is in the ez library
> # this produces the Greenhouse-Geisser corrected p-values
> 
> # We have to make sure that ezANOVA will recognize the group inputs as factors which we can do using recode
> # otherwise we get incorrect dfs and also would not obtain the Greenhouse-Geisser output
> 
> # factor usually works at constructing factors but doesn't seem to work with ezANOVA
> 
> #mdata$mood2 <- factor(mdata$mood, labels = c('A','B','C'))
> #mdata$time <- factor(mdata$time, labels=c('base','follow-up'))
> #mdata$group2 <- factor(mdata$group, labels=c('controls','patients'))
> 
> #  so we use recode which does
> 
> mdata$mood2 <- recode(mdata$mood, "1='A'; 2='B'; 3='C'")
> mdata$time <- recode(mdata$time, "1='base';2='follow-up'")
> #mdata$group2 <- recode(mdata$group, "1='controls';2='patients'")
> mdata$group2 <- recode(mdata$group2, "1='controls';2='patients'")
> 
> out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,type=3)
Warning: Converting "subject" to factor for ANOVA.
Warning: Converting "mood2" to factor for ANOVA.
Warning: Converting "time" to factor for ANOVA.
Warning: Converting "group2" to factor for ANOVA.
> out
$ANOVA
             Effect DFn DFd          SSn      SSd            F            p
1       (Intercept)   1  16 1020.5925926 18.03704 905.33059548 1.637243e-15
2            group2   1  16   13.3703704 18.03704  11.86036961 3.336341e-03
3             mood2   2  32   19.9074074 45.96296   6.92989525 3.158542e-03
5              time   1  16   23.1481481 37.37037   9.91080278 6.219390e-03
4      group2:mood2   2  32    0.1296296 45.96296   0.04512490 9.559388e-01
6       group2:time   1  16    0.1481481 37.37037   0.06342914 8.043623e-01
7        mood2:time   2  32    2.9074074 61.74074   0.75344931 4.789078e-01
8 group2:mood2:time   2  32    0.6851852 61.74074   0.17756449 8.381263e-01
  p<.05         ges
1     * 0.862202753
2     * 0.075760756
3     * 0.108772640
5     * 0.124279181
4       0.000794101
6       0.000907441
7       0.017512549
8       0.004183154

$`Mauchly's Test for Sphericity`
             Effect         W          p p<.05
3             mood2 0.9848657 0.89192368      
4      group2:mood2 0.9848657 0.89192368      
7        mood2:time 0.7082474 0.07522939      
8 group2:mood2:time 0.7082474 0.07522939      

$`Sphericity Corrections`
             Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF]
3             mood2 0.9850914 0.003325405         * 1.1223634 0.003158542
4      group2:mood2 0.9850914 0.954246131           1.1223634 0.955938799
7        mood2:time 0.7741421 0.449280837           0.8414513 0.458946950
8 group2:mood2:time 0.7741421 0.782942944           0.8414513 0.801527298
  p[HF]<.05
3         *
4          
7          
8          

install.packages("afex")
library(afex)
> 
# extracting residuals see https://stackoverflow.com/questions/30238213/ggplot2-residuals-with-ezanova

> anova_2 <- aov_ez("subject","score",mdata,within=c("mood2","time"),between=c("group2"))
Converting to factor: group2
Contrasts set to contr.sum for the following variables: group2
> nice(anova_2)
Anova Table (Type 3 tests)

Response: score
             Effect          df  MSE        F   ges p.value
1            group2       1, 16 1.13 11.86 **  .076    .003
2             mood2 1.97, 31.52 1.46  6.93 **  .109    .003
3      group2:mood2 1.97, 31.52 1.46     0.05 <.001    .954
4              time       1, 16 2.34  9.91 **  .124    .006
5       group2:time       1, 16 2.34     0.06 <.001    .804
6        mood2:time 1.55, 24.77 2.49     0.75  .018    .449
7 group2:mood2:time 1.55, 24.77 2.49     0.18  .004    .783
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Sphericity correction method: GG 
> residuals(anova_2)
Data was changed during ANOVA calculation. Thus, residuals cannot be added to original data.
residuals(..., append = TRUE) will return data and residuals.
         1          2          3          4          5          6          7 
 1.1111111 -1.3333333 -0.3333333  0.3333333 -0.6666667 -0.4444444  0.6666667 
         8          9         10         11         12         13         14 
-1.3333333 -1.8888889  0.6666667 -0.5555556  2.8888889 -0.3333333 -0.3333333 
        15         16         17         18         19         20         21 
-0.8888889 -0.3333333  0.4444444 -0.1111111 -1.3333333  0.6666667  0.1111111 
        22         23         24         25         26         27         28 
-1.3333333  1.4444444  0.8888889 -0.3333333  1.6666667  1.1111111  2.6666667 
        29         30         31         32         33         34         35 
 0.4444444 -1.1111111  0.6666667 -0.3333333 -0.8888889  1.6666667 -1.5555556 
        36         37         38         39         40         41         42 
-0.1111111 -1.3333333  0.6666667  0.1111111 -2.3333333  1.4444444 -2.1111111 
        43         44         45         46         47         48         49 
-0.3333333 -1.3333333  1.1111111 -1.3333333  0.4444444 -1.1111111  0.6666667 
        50         51         52         53         54         55         56 
-0.3333333  2.1111111 -0.3333333 -0.5555556 -0.1111111  1.6666667  0.6666667 
        57         58         59         60         61         62         63 
-0.8888889  0.6666667 -1.5555556  0.8888889  0.1111111 -0.3333333 -0.3333333 
        64         65         66         67         68         69         70 
-1.6666667 -1.6666667 -1.4444444 -0.8888889 -1.3333333 -0.3333333  0.3333333 
        71         72         73         74         75         76         77 
-0.6666667  1.5555556 -0.8888889 -0.3333333 -0.3333333  2.3333333 -0.6666667 
        78         79         80         81         82         83         84 
 0.5555556  0.1111111  2.6666667 -1.3333333  1.3333333  0.3333333 -0.4444444 
        85         86         87         88         89         90         91 
 1.1111111  1.6666667  0.6666667 -1.6666667  2.3333333 -1.4444444  0.1111111 
        92         93         94         95         96         97         98 
-0.3333333  1.6666667 -0.6666667  1.3333333 -2.4444444 -0.8888889  0.6666667 
        99        100        101        102        103        104        105 
-0.3333333 -2.6666667  0.3333333  1.5555556  0.1111111 -1.3333333  0.6666667 
       106        107        108 
 2.3333333 -0.6666667  2.5555556 
> 

>>>>>>>>>>>>>> ez.glm OBSOLETE

install.packages("afex")
library(afex)
out <- ezANOVA(data=mdata, dv=.(score), wid=.(subject), within=.(mood2,time),between=.(group2),detailed=T,type=3)
anova_1 <- ez.glm("subject","score",mdata,within=c("mood2","time"),between=c("group2"),return="lm")
> residuals(anova_1)
       A_base A_follow.up     B_base B_follow.up     C_base C_follow.up
1   1.1111111  -1.3333333 -0.3333333   0.3333333 -0.6666667  -0.4444444
2   0.1111111  -0.3333333 -0.3333333  -1.6666667 -1.6666667  -1.4444444
3  -0.8888889  -1.3333333 -0.3333333   0.3333333 -0.6666667   1.5555556
4  -0.8888889  -0.3333333 -0.3333333   2.3333333 -0.6666667   0.5555556
5   0.1111111   2.6666667 -1.3333333   1.3333333  0.3333333  -0.4444444
6   1.1111111   1.6666667  0.6666667  -1.6666667  2.3333333  -1.4444444
7   0.1111111  -0.3333333  1.6666667  -0.6666667  1.3333333  -2.4444444
8  -0.8888889   0.6666667 -0.3333333  -2.6666667  0.3333333   1.5555556
9   0.1111111  -1.3333333  0.6666667   2.3333333 -0.6666667   2.5555556
10  0.6666667  -1.3333333 -1.8888889   0.6666667 -0.5555556   2.8888889
11 -0.3333333  -0.3333333 -0.8888889  -0.3333333  0.4444444  -0.1111111
12 -1.3333333   0.6666667  0.1111111  -1.3333333  1.4444444   0.8888889
13 -0.3333333   1.6666667  1.1111111   2.6666667  0.4444444  -1.1111111
14  0.6666667  -0.3333333 -0.8888889   1.6666667 -1.5555556  -0.1111111
15 -1.3333333   0.6666667  0.1111111  -2.3333333  1.4444444  -2.1111111
16 -0.3333333  -1.3333333  1.1111111  -1.3333333  0.4444444  -1.1111111
17  0.6666667  -0.3333333  2.1111111  -0.3333333 -0.5555556  -0.1111111
18  1.6666667   0.6666667 -0.8888889   0.6666667 -1.5555556   0.8888889
> rstandard(anova_1)
       A_base A_follow.up     B_base B_follow.up C_base C_follow.up
1   1.3130643  -1.1547005 -0.3253957   0.2085144   -0.6 -0.30151134
2   0.1313064  -0.2886751 -0.3253957  -1.0425721   -1.5 -0.97991187
3  -1.0504515  -1.1547005 -0.3253957   0.2085144   -0.6  1.05528971
4  -1.0504515  -0.2886751 -0.3253957   1.4596009   -0.6  0.37688918
5   0.1313064   2.3094011 -1.3015827   0.8340577    0.3 -0.30151134
6   1.3130643   1.4433757  0.6507914  -1.0425721    2.1 -0.97991187
7   0.1313064  -0.2886751  1.6269784  -0.4170288    1.2 -1.65831240
8  -1.0504515   0.5773503 -0.3253957  -1.6681153    0.3  1.05528971
9   0.1313064  -1.1547005  0.6507914   1.4596009   -0.6  1.73369023
10  0.7878386  -1.1547005 -1.8439089   0.4170288   -0.5  1.95982374
11 -0.3939193  -0.2886751 -0.8677218  -0.2085144    0.4 -0.07537784
12 -1.5756772   0.5773503  0.1084652  -0.8340577    1.3  0.60302269
13 -0.3939193   1.4433757  1.0846523   1.6681153    0.4 -0.75377836
14  0.7878386  -0.2886751 -0.8677218   1.0425721   -1.4 -0.07537784
15 -1.5756772   0.5773503  0.1084652  -1.4596009    1.3 -1.43217889
16 -0.3939193  -1.1547005  1.0846523  -0.8340577    0.4 -0.75377836
17  0.7878386  -0.2886751  2.0608393  -0.2085144   -0.5 -0.07537784
18  1.9695965   0.5773503 -0.8677218   0.4170288   -1.4  0.60302269