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
```