Using aov()

In this section we will use a one-way, within-subjects ANOVA to evaluate the hypothesis that weight did not vary across days. First, we use aov. The random factor Chick is enclosed in Error, and we have also informed aov that the fixed factor, Day, is nested within Chick.

# check for balance!
# xtabs(~Chick+Day,ChickWeight) # check balance:
# Should we worry about sphericity? 
# Answer: yes, because correlations won't be constant for all days
chick.aov.01 <- aov(Weight~Day+Error(Chick/Day),data=ChickWeight)
summary(chick.aov.01) # assumes sphericity
## 
## Error: Chick
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19  12220   643.2               
## 
## Error: Chick:Day
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Day        5   2556   511.2   3.616 0.00488 **
## Residuals 95  13430   141.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using aov_car()

The \(p\) value for Day assumes that the variance-covariance matrix for the residuals is spherical. This condition is met when the correlations between the various levels of this Day are equally correlated. Note that this condition probably is not met in this case because we expect the correlation to be higher for days that are close together. So, we need to adjust the \(p\) value for possible violations of the sphericity assumpition. Here we use aov_car to perform the ANOVA, use the Mauchly test to evaluate sphericity, and adjust the \(p\) value.

library(afex) # afex library:
chick.aov.car.01  <- aov_car(Weight~Day+Error(Chick/Day),data=ChickWeight)
summary(chick.aov.car.01 )
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##             Sum Sq num Df Error SS den Df  F value    Pr(>F)    
## (Intercept) 634089      1    12220     19 985.8665 < 2.2e-16 ***
## Day           2556      5    13430     95   3.6164  0.004878 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##     Test statistic    p-value
## Day      0.0057318 1.2874e-12
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##      GG eps Pr(>F[GG])  
## Day 0.39602      0.037 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        HF eps Pr(>F[HF])
## Day 0.4418627 0.03159921

Using lmer()

Next, we analyze the data with a mixed model fit to the data with lmer in the lmerTest package. The random factor is specificed as (1 | Chick), which allows the model to adjust the intercept on a chick-by-chick basis. The model assumes that residuals for each level of the within-chick factor are independent, and therefore the \(F\) test for the fixed effect assumes sphericity. A statistical test on the random factor Chick is done with rand or ranova.

library(lmerTest) 
chick.lmer.01 <- lmer(Weight~ Day + (1 | Chick),data=ChickWeight)
anova(chick.lmer.01) # assumes sphericity
## Type III Analysis of Variance Table with Satterthwaite's method
##     Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Day 2556.1  511.23     5    95  3.6164 0.004878 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rand(chick.lmer.01) # same result as next line...
ranova(chick.lmer.01)  # Note that the p value is conservative
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## Weight ~ Day + (1 | Chick)
##             npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>         8 -469.16 954.32                         
## (1 | Chick)    7 -481.26 976.51 24.196  1    8.7e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The following block of code illustrates a method of evaluating the within-subject factor that does not rely on the assumption of sphericity. We first fit the data with a mixed model that does not include Day, and then use anova that compares the goodness-of-fit of the full model, which includes Day to the reduced model. The test is significant, so we reject the null hypothesis that the effects of Day are zero.

# alternative test of fixed effect:
# doesn't assume sphericity
chick.lmer.00 <- lmer(Weight~ 1 + (1 | Chick),data=ChickWeight)
# note that chi-square test tends to be conservative:
anova(chick.lmer.00,chick.lmer.01) # effect of Day is significant
## Data: ChickWeight
## Models:
## chick.lmer.00: Weight ~ 1 + (1 | Chick)
## chick.lmer.01: Weight ~ Day + (1 | Chick)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## chick.lmer.00    3 982.28 990.64 -488.14   976.28                        
## chick.lmer.01    8 974.85 997.15 -479.43   958.85 17.424  5   0.003763 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using lme()

Finally, we analyze the data with mixed models using lme in the nlme package. The advantage of the lme command is that it allows us to specify the form of the residuals variance-covariance matrix that we want to use in our model. The first model assumes that the pairwise correlations of the levels of Day are zero. The second model assumes that the correlations are non-zero but constant for all pairs of Days. In both cases, the the variance-covariance matrix is compound symmetric. Note that the value of 0.5 in the second lme call simply initializes the fitting procedure; the final estimated correlation, Rho, is 0.293.

library(nlme) 
chick.lme.01 <- lme(Weight~Day,
                    random=~1|Chick,
                    correlation=NULL,
                    data=ChickWeight)
chick.lme.02 <- lme(Weight~Day,
                    random=~1|Chick,
                    correlation=corCompSymm(value=0.5,form=~1|Chick),
                    data=ChickWeight)
anova(chick.lme.01) # assumes sphericity
##             numDF denDF  F-value p-value
## (Intercept)     1    95 985.8665  <.0001
## Day             5    95   3.6164  0.0049
anova(chick.lme.02) # assumes sphericity
##             numDF denDF  F-value p-value
## (Intercept)     1    95 985.8632  <.0001
## Day             5    95   3.6164  0.0049
summary(chick.lme.01) # [Rho, correlation between levels of Day, is zero and not shown]
## Linear mixed-effects model fit by REML
##   Data: ChickWeight 
##        AIC      BIC    logLik
##   954.3151 976.2047 -469.1576
## 
## Random effects:
##  Formula: ~1 | Chick
##         (Intercept) Residual
## StdDev:     9.14527 11.88966
## 
## Fixed effects:  Weight ~ Day 
##                Value Std.Error DF   t-value p-value
## (Intercept) 72.69167  2.315131 95 31.398511  0.0000
## Day1        -4.01855  2.426967 95 -1.655790  0.1011
## Day2        -4.22446  2.426967 95 -1.740634  0.0850
## Day3        -3.45538  2.426967 95 -1.423743  0.1578
## Day4        -0.41129  2.426967 95 -0.169467  0.8658
## Day5         4.00780  2.426967 95  1.651360  0.1020
##  Correlation: 
##      (Intr) Day1 Day2 Day3 Day4
## Day1  0.0                      
## Day2  0.0   -0.2               
## Day3  0.0   -0.2 -0.2          
## Day4  0.0   -0.2 -0.2 -0.2     
## Day5  0.0   -0.2 -0.2 -0.2 -0.2
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.02576628 -0.62023408 -0.05510695  0.49794771  2.09560207 
## 
## Number of Observations: 120
## Number of Groups: 20
summary(chick.lme.02) # note the value of Rho [correlation between levels of Day]
## Linear mixed-effects model fit by REML
##   Data: ChickWeight 
##        AIC      BIC    logLik
##   956.3151 980.9409 -469.1576
## 
## Random effects:
##  Formula: ~1 | Chick
##         (Intercept) Residual
## StdDev:    4.994437 14.14411
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | Chick 
##  Parameter estimate(s):
##       Rho 
## 0.2933777 
## Fixed effects:  Weight ~ Day 
##                Value Std.Error DF   t-value p-value
## (Intercept) 72.69167  2.315135 95 31.398458  0.0000
## Day1        -4.01855  2.426966 95 -1.655790  0.1011
## Day2        -4.22446  2.426966 95 -1.740635  0.0850
## Day3        -3.45538  2.426966 95 -1.423743  0.1578
## Day4        -0.41129  2.426966 95 -0.169467  0.8658
## Day5         4.00780  2.426966 95  1.651360  0.1020
##  Correlation: 
##      (Intr) Day1 Day2 Day3 Day4
## Day1  0.0                      
## Day2  0.0   -0.2               
## Day3  0.0   -0.2 -0.2          
## Day4  0.0   -0.2 -0.2 -0.2     
## Day5  0.0   -0.2 -0.2 -0.2 -0.2
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.08531590 -0.59078519 -0.03993286  0.53120711  2.14856901 
## 
## Number of Observations: 120
## Number of Groups: 20

Next, we fit a model that does not make any assumptions about the correlations between levels of Day. (Note that the variance-covariance matrix is symmetrical (correlation = corSymm) because the correlation between Days 1 and 2 is the same as the correlation between Days 2 and 1).

chick.lme.03 <- lme(Weight~Day,
                    random=~1|Chick,
                    correlation=corSymm(form=~1|Chick),
                    data=ChickWeight)
anova(chick.lme.03) # does NOT assume sphericity:
##             numDF denDF   F-value p-value
## (Intercept)     1    95 1416.8889  <.0001
## Day             5    95   11.0976  <.0001
summary(chick.lme.03) # note the values of Rho [correlation between levels of Day]
## Linear mixed-effects model fit by REML
##   Data: ChickWeight 
##        AIC      BIC    logLik
##   880.6963 943.6288 -417.3481
## 
## Random effects:
##  Formula: ~1 | Chick
##         (Intercept) Residual
## StdDev:    8.675637 12.23658
## 
## Correlation Structure: General
##  Formula: ~1 | Chick 
##  Parameter estimate(s):
##  Correlation: 
##   1      2      3      4      5     
## 2  0.205                            
## 3 -0.457  0.461                     
## 4 -0.625 -0.047  0.564              
## 5 -0.843 -0.360  0.411  0.775       
## 6 -0.845 -0.453  0.306  0.778  0.968
## Fixed effects:  Weight ~ Day 
##                Value Std.Error DF   t-value p-value
## (Intercept) 72.69167  2.315135 95 31.398453  0.0000
## Day1        -4.01855  3.603760 95 -1.115098  0.2676
## Day2        -4.22446  2.659087 95 -1.588689  0.1155
## Day3        -3.45538  1.838251 95 -1.879708  0.0632
## Day4        -0.41129  1.727191 95 -0.238127  0.8123
## Day5         4.00780  2.053052 95  1.952115  0.0539
##  Correlation: 
##      (Intr) Day1   Day2   Day3   Day4  
## Day1 -0.425                            
## Day2 -0.096  0.426                     
## Day3  0.295 -0.411  0.244              
## Day4  0.363 -0.672 -0.612 -0.026       
## Day5  0.176 -0.702 -0.831 -0.162  0.540
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.11744147 -0.77538428 -0.04953608  0.68534481  2.39658548 
## 
## Number of Observations: 120
## Number of Groups: 20
#anova(chick.lme.02,chick.lme.03) # model 2 is better

Finally, we can use anova to compare the fits attained by the three models. Note the difference in degrees of freedom. The first model estimates 8 coefficients: one intercept, 5 coefficients for the fixed effect Day, and the variance components for Chick and Residual. The second model estimates those 8 coefficients plus a ninth coefficient for the correlation (Rho) between levels of Day. The third model estimates 23 coefficients: an intercept, 5 effects for Day, variance components for Chick and Residuals, and 15 correlations between the levels of Day. Hence, the difference between the number of free parameters is \(23-9=14\). When the models are nested (as they are here), the test statistic is \(-2 \times (\mathit{logLik}_2 - \mathit{logLik}_1)\), or \(-2 \times -51.82 = 103.64\). When the null hypothesis of no difference between models is true, then the likelihood ratio statistic is distributed approximately as a chi-square variable with 14 degrees of freedom. In this case the third model provides a significantly better fit to the data.

anova(chick.lme.01,chick.lme.02,chick.lme.03) # model 3 is better
##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## chick.lme.01     1  8 954.3151 976.2047 -469.1576                        
## chick.lme.02     2  9 956.3151 980.9409 -469.1576 1 vs 2   0.0000       1
## chick.lme.03     3 23 880.6963 943.6288 -417.3481 2 vs 3 103.6188  <.0001

Variance Components & Intraclass Correlation for Chick.

The following code uses the ANOVA table to calculate the ANOVA estimates of the variance components and intraclass correlations for Chick and Residuals.

summary(chick.aov.01)
## 
## Error: Chick
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19  12220   643.2               
## 
## Error: Chick:Day
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Day        5   2556   511.2   3.616 0.00488 **
## Residuals 95  13430   141.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
var.comp.resid <- 141.36
( var.comp.chick <- (643.2 - var.comp.resid) / 6 )# denominator equals number of levels of Day
## [1] 83.64
( resid.icc <- var.comp.resid / (var.comp.chick + var.comp.resid) )
## [1] 0.6282667
( chick.icc <- var.comp.chick / (var.comp.chick + var.comp.resid) )
## [1] 0.3717333

These values also can be calculated with anovaMM in the VCA package:

# ANOVA estimates of var comp:
library(VCA)
aov.vca <- anovaMM(Weight~Day+(Chick),Data=ChickWeight)
print(aov.vca,digits=3)
## 
## 
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
## 
##  [Fixed Effects]
## 
##    int  Dayd0 Dayd10  Dayd2  Dayd4  Dayd6  Dayd8 
## 76.699 -8.026  4.094 -8.232 -7.463 -4.419  0.000 
## 
## 
##  [Variance Components]
## 
##   Name  DF     SS        MS      VC      %Total SD    CV[%] 
## 1 total 67.421                   225     100    15    20.635
## 2 Chick 19     12220.416 643.18  83.636  37.172 9.145 12.581
## 3 error 95     13429.584 141.364 141.364 62.828 11.89 16.356
## 
## Mean: 72.692 (N = 120) 
## 
## Experimental Design: balanced  |  Method: ANOVA

The variance components estimated from the model fit by lmer are equal to the anova estimates.

# using lmer
chick.vca <- VarCorr(chick.lmer.01)
print(chick.vca,comp=c("Variance","Std.Dev."))
##  Groups   Name        Variance Std.Dev.
##  Chick    (Intercept)  83.636   9.1453 
##  Residual             141.364  11.8897
library(performance)
icc(chick.lmer.01,by_group=T)
## # ICC by Group
## 
## Group |   ICC
## -------------
## Chick | 0.372

The values of the variance components depend on the structure of the residuals variance-covariance matrix. Recall that the third model fit the data best, and therefore we might favor the variance components estimated by that model over the other estimates.

# using lme:
VarCorr(chick.lme.01) # model assumes compound symmetry [rho = 0]
## Chick = pdLogChol(1) 
##             Variance  StdDev  
## (Intercept)  83.63596  9.14527
## Residual    141.36404 11.88966
VarCorr(chick.lme.02) # model assumes compound symmetry [rho = 0.293]
## Chick = pdLogChol(1) 
##             Variance StdDev   
## (Intercept)  24.9444  4.994437
## Residual    200.0559 14.144111
VarCorr(chick.lme.03) # model does NOT assume compound symmetry
## Chick = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept)  75.26668  8.675637
## Residual    149.73386 12.236579

Association Strength & Effect Size for Day.

Measures of association strength and effect size for Day are calculated the normal way with commands in the effectsize package.

# effect size
library(effectsize)
cohens_f(chick.aov.01)
## # Effect Size for ANOVA (Type I)
## 
## Group     | Parameter | Cohen's f (partial) |      95% CI
## ---------------------------------------------------------
## Chick:Day |       Day |                0.44 | [0.18, Inf]
## 
## - One-sided CIs: upper bound fixed at [Inf].
cohens_f(chick.aov.car.01 )
## # Effect Size for ANOVA (Type III)
## 
## Parameter | Cohen's f (partial) |      95% CI
## ---------------------------------------------
## Day       |                0.44 | [0.18, Inf]
## 
## - One-sided CIs: upper bound fixed at [Inf].
cohens_f(chick.lmer.01)
## # Effect Size for ANOVA (Type III)
## 
## Parameter | Cohen's f (partial) |      95% CI
## ---------------------------------------------
## Day       |                0.44 | [0.18, Inf]
## 
## - One-sided CIs: upper bound fixed at [Inf].

# association strength
library(effectsize)
omega_squared(chick.aov.01)
## # Effect Size for ANOVA (Type I)
## 
## Group     | Parameter | Omega2 (partial) |       95% CI
## -------------------------------------------------------
## Chick:Day |       Day |             0.06 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(chick.aov.car.01 )
## # Effect Size for ANOVA (Type III)
## 
## Parameter | Omega2 (partial) |       95% CI
## -------------------------------------------
## Day       |             0.06 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].
omega_squared(chick.lmer.01)
## # Effect Size for ANOVA (Type III)
## 
## Parameter | Omega2 (partial) |       95% CI
## -------------------------------------------
## Day       |             0.11 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Linear Contrasts

In this section we illustrate how to evaluate linear contrasts on the within-chick variable.

The function contr.poly can be used to create weights for polynomial contrasts – i.e., the contrast weights needed for analyses of trends. First we use contr.poly to create weights that can be used to test for linear and quadratic trends of weight across days and verify that the two sets of weights are orthogonal.

# levels of day are space equally, so we can use default contr.poly:
contr.poly(n=6)
##              .L         .Q         .C         ^4          ^5
## [1,] -0.5976143  0.5455447 -0.3726780  0.1889822 -0.06299408
## [2,] -0.3585686 -0.1091089  0.5217492 -0.5669467  0.31497039
## [3,] -0.1195229 -0.4364358  0.2981424  0.3779645 -0.62994079
## [4,]  0.1195229 -0.4364358 -0.2981424  0.3779645  0.62994079
## [5,]  0.3585686 -0.1091089 -0.5217492 -0.5669467 -0.31497039
## [6,]  0.5976143  0.5455447  0.3726780  0.1889822  0.06299408
lin.weights <- contr.poly(n=6)[,1] # 1st column
quad.weights <- contr.poly(n=6)[,2] # second column
sum(lin.weights*quad.weights) # orthogonal!
## [1] -5.551115e-17

Now we use the linear weights to convert the data from each chick to a single variable that represents the linear trend score, and then use a \(t\) test to evaluate the null hypothesis that the linear trend of Weight across Day is zero. Note that we use the wide-format data frame cw.wide in the code below. Note that the degrees of freedom equals \(N-1=19\).

head(cw.wide) # data in wide format
##   Chick       d0       d2       d4       d6       d8      d10
## 1    c1 83.65735 80.42716 53.00276 64.37642 44.90555 50.50263
## 2    c2 83.65735 57.64629 58.41394 58.41117 65.70019 70.78673
## 3    c3 71.66997 97.51282 85.46982 73.32430 71.17246 68.75832
## 4    c4 95.64474 80.42716 69.23629 61.39379 71.17246 74.16741
## 5    c5 47.69519 46.25585 53.00276 61.39379 67.88910 72.13900
## 6    c6 47.69519 46.25585 53.00276 61.39379 73.36137 81.60491
y.mat <- as.matrix(cw.wide[,-1]) # extract all columns except 1st one and save as matrix
lin.scores <- y.mat %*% lin.weights; # linear trend scores for each chick
t.test(lin.scores)
## 
##  One Sample t-test
## 
## data:  lin.scores
## t = 2.2646, df = 19, p-value = 0.03543
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   0.7999345 20.3180804
## sample estimates:
## mean of x 
##  10.55901

Here we use aov to evaluate the linear and non-linear trends. Note that the degrees of freedom equals \((a-1)(n-1)=5\times 19=95\).

# use aov/split:
ChickWeight$Day <- as.ordered(ChickWeight$Day)
class(ChickWeight$Day) # an >>ordered<< factor
## [1] "ordered" "factor"
contrasts(ChickWeight$Day) # default contrast are POLYNOMIAL trends (linear, quad, etc)
##              .L         .Q         .C         ^4          ^5
## [1,] -0.5976143  0.5455447 -0.3726780  0.1889822 -0.06299408
## [2,] -0.3585686 -0.1091089  0.5217492 -0.5669467  0.31497039
## [3,] -0.1195229 -0.4364358  0.2981424  0.3779645 -0.62994079
## [4,]  0.1195229 -0.4364358 -0.2981424  0.3779645  0.62994079
## [5,]  0.3585686 -0.1091089 -0.5217492 -0.5669467 -0.31497039
## [6,]  0.5976143  0.5455447  0.3726780  0.1889822  0.06299408
chick.aov.01b <- aov(Weight~Day+Error(Chick/Day),data=ChickWeight)
# linear trend is significant; nonlinear trends are not
summary(chick.aov.01b,split=list(Day=list(linear=1,others=2:5)))
## 
## Error: Chick
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19  12220   643.2               
## 
## Error: Chick:Day
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Day            5   2556   511.2   3.616 0.004878 ** 
##   Day: linear  1   2230  2229.9  15.774 0.000139 ***
##   Day: others  4    326    81.6   0.577 0.679986    
## Residuals     95  13430   141.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
F.linear <- 15.77
sqrt(F.linear) # convert to t statistic to compare to next test
## [1] 3.971146

The next two chunks of code use emmeans to analyze the linear trend. The first chunk uses an aov object. Note that the degrees of freedom is 95.

# emmeans uses GLOBAL error (df=95) here:
library(emmeans)
chick.emm <- emmeans(chick.aov.01,specs=~Day)
contrast(chick.emm,method=list(linear=lin.weights))
##  contrast estimate   SE df t.ratio p.value
##  linear       10.6 2.66 95   3.972  0.0001

Now we do the same analysis on an aov_car object. Here the degrees of freedom equals 19, which is the same value that we had in our \(t\) test.

# emmeans uses LOCAL error (df=19) here:
options(contrasts=c("contr.sum","contr.poly"))
# use aov_car !
chick.aov.01c <- aov_car(Weight~Day+Error(Chick/Day),data=ChickWeight)
library(emmeans)
chick.aov.car.emm <- emmeans(chick.aov.01c,specs=~Day)
contrast(chick.aov.car.emm,method=list(linear=lin.weights))
##  contrast estimate   SE df t.ratio p.value
##  linear       10.6 4.66 19   2.265  0.0354

Questions (trackbox data)

The data for this section can be loaded into R with the following command:

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/trackboxData.rda"))

Discovery Day is a day set aside by the United States Naval Postgraduate School in Monterey, California, to invite the general public into its laboratories. On Discovery Day, 21 October 1995, data on reaction time and hand-eye coordination were collected on 118 members of the public who visited the Human Systems Integration Laboratory. The age and sex of each subject were also recorded. Visitors were mostly in family groups. A rotary pursuit tracking experiment was done to examine motor learning and hand-eye coordination. The equipment was a rotating disk with a 3/4” target spot. In the Circle condition, the target spot moved at a constant speed in a circular path. In the Box condition, the target spot moved at various speeds as it moved along a box-shaped path. The subject’s task was to maintain contact with the target spot with a metal wand. Four trials were recorded for each of 108 subjects. Each trial lasted 15~s, and the total contact duration during each trial was recorded. The data from the Box condition (n=70) are stored in a wide-format data frame, track.wide. The variables are sex, age, subject, time1, time2, time3, time4. The dependent variable, the amount of time the subject maintained contact with the target on each trial, is in the variables time1-4. The data also are stored in the long-format data frame, track.long, which contains the variables sex, age, subject, times, and contact. Ignore the factors sex and age in the following analyses.

1. Use a within-subjects (repeated-measures) ANOVA to determine if contact time differed across trials. Adjust the \(p\) value for the effect of trial with the Huynh-Feldt correction for sphericity.

2. Calculate the effect size of times.

3. Calculate the variance component and intraclass correlation for subjects.

4. Evaluate the null hypothesis that the variance component for subjects is zero.

5. Use a linear contrast to test for a positive (i.e., increasing) linear trend in contact time across trials.