Protected Trend Analysis

Tensile Strength Example

> # The strength were coded. Add 15 to get original data.

> # rep(5,5) = repeat 5 five times (2nd argument=repetitions)

> # seq(15,35,5) = sequence from 15 to 35 with increment 5

> # rep(seq(15,35,5), rep(5,5)) = 15 five times, 20 five times,

> # ..., 35 five times

> fiber <- data.frame(strength=scan()+15,

+ cotton=rep( seq(15,35,5), rep(5,5) ) )

1: -8 -8 0 -4 -6 -3 2 -3 3 3 -1 3 3 4 4

16: 4 10 7 4 8 -8 -5 -4 0 -4

26:

Read 25 items

> ybars <- by(fiber$strength, fiber$cotton, mean)

> ybars # treatment means

INDICES: 15
[1] 9.8
--------------------------------------------------------- 
INDICES: 20
[1] 15.4
--------------------------------------------------------- 
INDICES: 25
[1] 17.6
--------------------------------------------------------- 
INDICES: 30
[1] 21.6
--------------------------------------------------------- 
INDICES: 35
[1] 10.8

> vars <- by(fiber$strength, fiber$cotton, var)

> vars # treatment variances

INDICES: 15
[1] 11.2
--------------------------------------------------------- 
INDICES: 20
[1] 9.8
--------------------------------------------------------- 
INDICES: 25
[1] 4.3
--------------------------------------------------------- 
INDICES: 30
[1] 6.8
--------------------------------------------------------- 
INDICES: 35
[1] 8.2

> # max variance = 11.2 < 4 × min variance

> # variances seem to be homogeneous, see Levene's test below

>

> # Store the folder for convenience.

> s664<-"C:\\Documents and Settings\\WMUFaculty\\My Documents\\s664\\"

> pdf(file=paste(s664,"fiber01.pdf",sep=""),

+ width=4,height=2.5,pointsize=10, version="1.4")

> par(mar=c(4.1,4.1,.5,.1), las=1)

> # Produces treatment means plot (main effect plot)

> plot(seq(15,35,5), as.vector(ybars), xlab="%Cotton",

+ ylab="mean tensile strength",pch=21,bg=rgb(0,.5,.05,.2))

> dev.off()

null device 
          1

> # Click fiber01.pdf to see graph

>

> # Levene's test of homoscedasticity

> car::leveneTest(strength ~ cotton, data=fiber)

Error in leveneTest.formula(strength ~ cotton, data=fiber) : 
 Levene's test is not appropriate with quantitative explanatory variables.

> is.factor(fiber[[2]])

[1] FALSE

> car::leveneTest(strength ~ factor(cotton), data=fiber)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  4  0.3179 0.8626
      20

> # The insignificant result indicates that the variances are homogeneous.

>

> # Now perform stage 1 F test of the Protected Trend Analysis

> summary(aov(strength~factor(cotton),data=fiber))

               Df Sum Sq Mean Sq F value    Pr(>F)    
factor(cotton)  4 475.76  118.94  14.757 9.128e-06 ***
Residuals      20 161.20    8.06                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

>

> # Significant result → Do stage 2 trend component tests

>

>

> # Using the distinct values of variable cotton,

> # poly(cotton,4) generates 4 orthogonal polynomial contrasts

> # Note: contr.poly(5) were used for this.

> print(fiber.lm.op<-lm(strength~poly(cotton,4),data=fiber))

Call:
lm(formula = strength ~ poly(cotton, 4), data = fiber)

Coefficients:
     (Intercept)  poly(cotton, 4)1  poly(cotton, 4)2  poly(cotton, 4)3  
          15.040             5.798           -18.526            -8.061  
poly(cotton, 4)4  
          -5.826  

> # Intercept = grand mean. The rest are contrasts' estimates, each

> # = square root of nLi × integer coeff contrast (as in notes).

> # (contrast)2 = Component

>

> # So, to get the estimates of the coefficients as in the class notes:

> x <- coef(fiber.lm.op)

> x[-1] <- x[-1] / sqrt(5*c(10,14,10,70))

> x # display β.hat's

     (Intercept) poly(cotton, 4)1 poly(cotton, 4)2 poly(cotton, 4)3
      15.0400000        0.8200000       -2.2142857       -1.1400000
poly(cotton, 4)4
      -0.3114286

> #

> summary(fiber.lm.op)

Call:
lm(formula = strength ~ poly(cotton, 4), data = fiber)

Residuals:
   Min     1Q Median     3Q    Max 
  -3.8   -2.6    0.4    1.4    5.2 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       15.0400     0.5678  26.488  < 2e-16 ***
poly(cotton, 4)1   5.7983     2.8390   2.042   0.0545 .  
poly(cotton, 4)2 -18.5260     2.8390  -6.526 2.33e-06 ***
poly(cotton, 4)3  -8.0610     2.8390  -2.839   0.0101 *  
poly(cotton, 4)4  -5.8263     2.8390  -2.052   0.0535 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 2.839 on 20 degrees of freedom
Multiple R-Squared: 0.7469,     Adjusted R-squared: 0.6963 
F-statistic: 14.76 on 4 and 20 DF,  p-value: 9.128e-06

>

> # Command below extracts coefficient (can use coef for short).

> coefficients(fiber.lm.op)

     (Intercept) poly(cotton, 4)1 poly(cotton, 4)2 poly(cotton, 4)3 
       15.040000         5.798276       -18.526043        -8.061017 
poly(cotton, 4)4 
       -5.826295

>

> msq <- coef(fiber.lm.op)[-1]^2

> msq # SS = MS = Components

poly(cotton, 4)1 poly(cotton, 4)2 poly(cotton, 4)3 poly(cotton, 4)4 
        33.62000        343.21429         64.98000         33.94571

>

> # Square the above t ratios for contrasts → F ratios

>

> # Can also do the following to get F ratios

> anova(fiber.lm.op)

Analysis of Variance Table

Response: strength
                Df Sum Sq Mean Sq F value    Pr(>F)    
poly(cotton, 4)  4 475.76  118.94  14.757 9.128e-06 ***
Residuals       20 161.20    8.06                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

>

> # Extract MSE from ANOVA table and find F ratios

> f0 <- msq / anova(fiber.lm.op)["Residuals", "Mean Sq"]

> f0 # F ratios of Components

poly(cotton, 4)1 poly(cotton, 4)2 poly(cotton, 4)3 poly(cotton, 4)4 
        4.171216        42.582418         8.062035         4.211627 

> pf(f0, 1, 20, lower=FALSE) # P-values of Components

poly(cotton, 4)1 poly(cotton, 4)2 poly(cotton, 4)3 poly(cotton, 4)4 
    5.452366e-02     2.325546e-06     1.013339e-02     5.346882e-02

> # P values are the same as those in t tests

>

> # The analysis above (see P values) indicates that it seems

> # to be appropriate to fit 3rd order polynomial regression.

>

> # poly(cotton, degree=3, raw=T) generates 3rd order polynomial

> # of cotton. The model formula is equivalent to

> # strength ~ cotton + I(cotton^2) + I(cotton^3)

> # where I is the 'as is' function

> summary(lm(strength~poly(cotton, degree=3, raw=T), data=fiber))

Call:
lm(formula = strength ~ poly(cotton, degree = 3, raw = T), data = fiber)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4686 -1.4686 -0.4686  2.6457  4.8886 

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)  
(Intercept)                        62.611429  39.757436   1.575   0.1302  
poly(cotton, degree = 3, raw = T)1 -9.011429   5.196609  -1.734   0.0976 .
poly(cotton, degree = 3, raw = T)2  0.481429   0.216046   2.228   0.0369 *
poly(cotton, degree = 3, raw = T)3 -0.007600   0.002874  -2.644   0.0152 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 3.048 on 21 degrees of freedom
Multiple R-Squared: 0.6936,     Adjusted R-squared: 0.6499 
F-statistic: 15.85 on 3 and 21 DF,  p-value: 1.295e-05