# We started by demonstrating that Helmert contrasts can
# shed light on Eysenck's theory of memory and level of processing.
# (See the Powerpoint for an explanation of Helmert contrasts.)
# First we implement these using a contrast coding system. The
# first, h1, compares the Intentional group with all of the others.
# The second, h2, compares Imagery (which should be the highest
# group) with the three groups that should be lower. The third (h3)
# compares the Adjective group with the two that should be lower, and
# finally h4 compares the two lowest levels of processing.

> h1 <- c(rep(-1,40),rep(4,10))
> h2 <- c(rep(-1,30),rep(3,10),rep(0,10))
> h3 <- c(rep(-1,20),rep(2,10),rep(0,20))
> h4 <- c(rep(-1,10),rep(1,10),rep(0,30))

# Note that we have a system of 4 variables that identify the
# five groups with unique patterns of values:

> cbind(h1,h2,h3,h4)
      h1 h2 h3 h4
 [1,] -1 -1 -1 -1  # Counting
 [2,] -1 -1 -1 -1
 [3,] -1 -1 -1 -1
 [4,] -1 -1 -1 -1
 [5,] -1 -1 -1 -1
 [6,] -1 -1 -1 -1
 [7,] -1 -1 -1 -1
 [8,] -1 -1 -1 -1
 [9,] -1 -1 -1 -1
[10,] -1 -1 -1 -1
[11,] -1 -1 -1  1  # Rhyming
[12,] -1 -1 -1  1
[13,] -1 -1 -1  1
[14,] -1 -1 -1  1
[15,] -1 -1 -1  1
[16,] -1 -1 -1  1
[17,] -1 -1 -1  1
[18,] -1 -1 -1  1
[19,] -1 -1 -1  1
[20,] -1 -1 -1  1
[21,] -1 -1  2  0  # Adjective
[22,] -1 -1  2  0
[23,] -1 -1  2  0
[24,] -1 -1  2  0
[25,] -1 -1  2  0
[26,] -1 -1  2  0
[27,] -1 -1  2  0
[28,] -1 -1  2  0
[29,] -1 -1  2  0
[30,] -1 -1  2  0
[31,] -1  3  0  0  # Imagery
[32,] -1  3  0  0
[33,] -1  3  0  0
[34,] -1  3  0  0
[35,] -1  3  0  0
[36,] -1  3  0  0
[37,] -1  3  0  0
[38,] -1  3  0  0
[39,] -1  3  0  0
[40,] -1  3  0  0
[41,]  4  0  0  0  # Intentional
[42,]  4  0  0  0
[43,]  4  0  0  0
[44,]  4  0  0  0
[45,]  4  0  0  0
[46,]  4  0  0  0
[47,]  4  0  0  0
[48,]  4  0  0  0
[49,]  4  0  0  0
[50,]  4  0  0  0
> attach(Eysenck)

# We get the same ANOVA, and, because we have an
# orthogonal contrast coding system, the t statistics
# test the four contrasts:

> summary(lm(Score~h1+h2+h3+h4))

Call:
lm(formula = Score ~ h1 + h2 + h3 + h4)

Residuals:
   Min     1Q Median     3Q    Max 
 -7.00  -1.85  -0.45   2.00   9.60 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.0600     0.4399  22.872  < 2e-16 ***
h1            0.4850     0.2199   2.205  0.03258 *     # Intentional learning differs from the other 4 groups.
h2            1.2750     0.2839   4.491 4.91e-05 ***   # Imagery differs from the average of the 3 lower groups.
h3            1.3500     0.4015   3.362  0.00159 **    # Adjective differs from the average of the 2 lower groups.
h4           -0.0500     0.6955  -0.072  0.94300       # The two lowest groups are not significantly different.
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.11 on 45 degrees of freedom
Multiple R-squared:  0.4468,    Adjusted R-squared:  0.3976 
F-statistic: 9.085 on 4 and 45 DF,  p-value: 1.815e-05


# We can also implement these contrasts by attaching the correct
# contrast matrix to the Group factor:

> #            Count  Adject  Image  Inten  Rhyme
> helmerts <- matrix(c(
+                -1,     -1,    -1,    4,      -1,
+                -1,     -1,     3,    0,      -1,
+                 2,     -1,     0,    0,      -1,
+                 0,     -1,     0,    0,       1), 4,5,byrow=TRUE)
> contrasts(Group) <- t(helmerts)

# This gives exactly the same result as our contrast
# coding system:

> summary(lm(Score~Group))

Call:
lm(formula = Score ~ Group)

Residuals:
   Min     1Q Median     3Q    Max 
 -7.00  -1.85  -0.45   2.00   9.60 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.0600     0.4399  22.872  < 2e-16 ***
Group1        0.4850     0.2199   2.205  0.03258 *  
Group2        1.2750     0.2839   4.491 4.91e-05 ***
Group3        1.3500     0.4015   3.362  0.00159 ** 
Group4       -0.0500     0.6955  -0.072  0.94300    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.11 on 45 degrees of freedom
Multiple R-squared:  0.4468,    Adjusted R-squared:  0.3976 
F-statistic: 9.085 on 4 and 45 DF,  p-value: 1.815e-05

# We can also create a Helmert contrast system by looking in
# the other direction: compare the lowest group with the three
# that should be above it; compare the next to lowest group with
# the two that should be above it (and so on).

> helmert2 <- matrix(c(
+                -1,     -1,    -1,     4,     -1,
+                 1,     -3,     1,     0,      1,
+                 1,      0,     1,     0,     -2,
+                 1,      0,    -1,     0,      0),4,5,byrow=TRUE)
> contrasts(Group) <- t(helmert2)
> contrasts(Group)
          [,1] [,2] [,3] [,4]
adjective   -1    1    1    1
counting    -1   -3    0    0
imagery     -1    1    1   -1
intent       4    0    0    0
rhyming     -1    1   -2    0

# That approach also provides support for Eysenck's theory:

> summary(lm(Score~Group))

Call:
lm(formula = Score ~ Group)

Residuals:
   Min     1Q Median     3Q    Max 
 -7.00  -1.85  -0.45   2.00   9.60 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.0600     0.4398  22.872  < 2e-16 ***
Group1        0.4850     0.2199   2.205  0.03258 *  
Group2        0.8583     0.2839   3.023  0.00412 ** 
Group3        1.7667     0.4015   4.400 6.58e-05 ***
Group4       -1.2000     0.6955  -1.725  0.09130 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.11 on 45 degrees of freedom
Multiple R-squared:  0.4468,    Adjusted R-squared:  0.3976 
F-statistic: 9.085 on 4 and 45 DF,  p-value: 1.815e-05

# We turned our attention to the subject of power. The
# simplest example possible would be a Z test about a
# single mean. If we work with an alpha level of .05, the
# critical value of the Z statistic is +-1.96. Here's
# what the sampling distribution of the statistic should
# look like if the null hypothesis is TRUE:

> x <- seq(-5,5,.01)
> y <- dnorm(x)
> plot(x,y,type='l',xlab="Z Statistic",ylab="",axes=FALSE)
> abline(h=0)
> axis(side=1,seq(-5,5,1))

# And here is the rejection region:

> lines(c(-1.96,-1.96),c(0,dnorm(-1.96)))
> lines(c(1.96,1.96),c(0,dnorm(1.96)))

# If the true mean is 110 instead of 100 (see Powerpoint),
# then we are actually sampling from THIS distribution:

> y2 <- dnorm(x, 10/3, 1)
> lines(x,y2)

# We'll still reject if we get a value below -1.96 (which
# is very unlikely)...

> pnorm(-1.96, 10/3, 1)
[1] 6.005335e-08

# ...or above 1.96:

> 1-pnorm(1.96, 10/3, 1)
[1] 0.9151756

# So the power is the sum of those two probabilities:

> 1-pnorm(1.96, 10/3, 1) +  pnorm(-1.96, 10/3, 1)
[1] 0.9151757

# More realistically, we would estimate sigma and use a t
# test. Here's the appropriate t distribution if the null
# hypothesis is TRUE:

> x <- seq(-8,8,.01)
> y <- dt(x, 24)
> plot(x,y,type='l',xlab="t statistic",ylab="",axes=FALSE); abline(h=0)
> axis(side=1,at=seq(-8,8,1))

# If the null hypothesis is false in the way we stipulated (mu=110),
# then we would be sampling from a noncentral t distribution with
# 24 degrees of freedom and noncentrality parameter = (110-100)/3 = 10/3:

> delta <- (10/3)

# Here's that distribution. Note the slight positive skew:

> y2 <- dt(x, 24, delta)
> lines(x,y2)

# We will reject the null whenever we get a test statistic that's
# bigger than +-2.06:

> tcrit <- qt(.975, 24)
> tcrit
[1] 2.063899

# So our power is:

> pt(-tcrit, 24, delta) + (1-pt(tcrit, 24, delta))
[1] 0.892017
>