Yage Ding

Rensselaer Polytechnic Institute

December 6, 2016

1. Setting

System under test

The Ecdat::Benefits data set resides in the Ecdat package of the R programming language. It contains a information of unemployement of blue collar workers.

library(Ecdat)
B <- Ecdat::Benefits

Factors and Levels

There are 18 parameter in this data set:

colnames(B)
##  [1] "stateur"  "statemb"  "state"    "age"      "tenure"   "joblost" 
##  [7] "nwhite"   "school12" "sex"      "bluecol"  "smsa"     "married" 
## [13] "dkids"    "dykids"   "yrdispl"  "rr"       "head"     "ui"

In this experiment, we choose factors sex, nwhite (color), age, and tenure (years of tenure in job lost) as independent variables, and stateur (state unemployment rate (in %)) to be the response variable. In the original data frame, sex and nwhite are 2-level factors, whileage and tenure are continuous variables.

B <- Ecdat::Benefits [c("sex", "nwhite", "age", "tenure", "stateur")]
head(B) 
##      sex nwhite age tenure stateur
## 1   male     no  49     21     4.5
## 2   male     no  26      2    10.5
## 3 female     no  40     19     7.2
## 4 female    yes  51     17     5.8
## 5   male     no  33      1     6.5
## 6   male     no  51      3     7.5
str(B)
## 'data.frame':    4877 obs. of  5 variables:
##  $ sex    : Factor w/ 2 levels "female","male": 2 2 1 1 2 2 2 2 2 2 ...
##  $ nwhite : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
##  $ age    : int  49 26 40 51 33 51 30 26 54 31 ...
##  $ tenure : int  21 2 19 17 1 3 5 3 20 1 ...
##  $ stateur: num  4.5 10.5 7.2 5.8 6.5 7.5 5.8 5.8 7.7 6 ...

We first transform age and tenure into 2 3-level factors.

##transform age into 3-level factor
for (i in 1:4877) {
  if (B$age[i] <= 30) {
  B$age[i] <- "<=30"
  } else if (B$age[i] > 40) {
    B$age[i] <- ">40"
    } else {
       B$age[i] <- ">30 & <=40"
    }
}

We then transform age, and tenure into factors. Now, we have a data frame with 4 categorical factors, and a continuous response variable. Factor sex and nwhite both have 2 levels, while age and tenure have 3.

B1 <- B[c("sex", "nwhite", "age", "tenure","stateur")]
B1$sex <- B$sex
B1$nwhite <- B$nwhite
B1$age <- factor(B$age)
B1$tenure <- factor(B$tenure)
str(B1)
## 'data.frame':    4877 obs. of  5 variables:
##  $ sex    : Factor w/ 2 levels "female","male": 2 2 1 1 2 2 2 2 2 2 ...
##  $ nwhite : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
##  $ age    : Factor w/ 3 levels "<=30",">30 & <=40",..: 3 1 2 3 2 3 1 1 3 2 ...
##  $ tenure : Factor w/ 3 levels "<=2.5",">2.5 & <=6",..: 3 1 1 1 1 2 2 2 2 1 ...
##  $ stateur: num  4.5 10.5 7.2 5.8 6.5 7.5 5.8 5.8 7.7 6 ...
##sex
levels(B1$sex)
## [1] "female" "male"
##nwhite
levels(B1$nwhite)
## [1] "no"  "yes"
##age
levels(B1$age)
## [1] "<=30"       ">30 & <=40" ">40"
##tenure
levels(B1$tenure)
## [1] "<=2.5"      ">2.5 & <=6" ">6"

Since we want to construct a \(2^k\) design, where all factors are limited to have only 2 levels. We converted each 3-level factor (age and tenure) to 2 2-level factors.

##Transform 3-level factor (age) into 2 2-level factors (age1 and age2)
for (i in 1:4877) {
  if (B$age[i]=="<=30"){
    B$age1[i]<-"30"
    B$age2[i]<-"30"
    } else if (B$age[i]==">40") {
      B$age2[i]<-"40"
      B$age1[i]<-"40"
      } else {
        a <- c(1,2)
        b <- c("30", "40")
        c<- sample(a,1)
        B$age1[i]<- b[c]
        B$age2[i]<- b[3-c]
        }
}

And now, our experiment has 6 2-level factors instead of 4 mixed level factors:

str(B2)
## 'data.frame':    4877 obs. of  7 variables:
##  $ sex    : Factor w/ 2 levels "female","male": 2 2 1 1 2 2 2 2 2 2 ...
##  $ nwhite : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
##  $ age1   : Factor w/ 2 levels "30","40": 2 1 2 2 1 2 1 1 2 2 ...
##  $ age2   : Factor w/ 2 levels "30","40": 2 1 1 2 2 2 1 1 2 1 ...
##  $ tenure1: Factor w/ 2 levels "2.5","6": 2 1 1 1 1 2 1 1 1 1 ...
##  $ tenure2: Factor w/ 2 levels "2.5","6": 2 1 1 1 1 1 2 2 2 1 ...
##  $ stateur: num  4.5 10.5 7.2 5.8 6.5 7.5 5.8 5.8 7.7 6 ...

To have a simpler and more clear view of the levels, we code the levels so that the lower levels are all represented as “-1”, while the higher levels are represented as “1”.

B$sex <- as.character(B$sex)
for (i in 1:4877){
  if (as.character(B$sex[i]) == "female") {
    B$sex[i] <- "-1"
  } else {
    B$sex[i] <- "1"
  }
}

Now, our design is a \(2^2 * 2^2 * 2^2 = 2^6\) deisgn, where there are 6 factors and each with 2 levels (-1 and 1). The final structure of our data frame is shown here:

str(B)
## 'data.frame':    4877 obs. of  7 variables:
##  $ sex    : Factor w/ 2 levels "-1","1": 2 2 1 1 2 2 2 2 2 2 ...
##  $ nwhite : Factor w/ 2 levels "-1","1": 1 1 1 2 1 1 1 1 1 1 ...
##  $ age1   : Factor w/ 2 levels "-1","1": 2 1 2 2 1 2 1 1 2 2 ...
##  $ age2   : Factor w/ 2 levels "-1","1": 2 1 1 2 2 2 1 1 2 1 ...
##  $ tenure1: Factor w/ 2 levels "-1","1": 2 1 1 1 1 2 1 1 1 1 ...
##  $ tenure2: Factor w/ 2 levels "-1","1": 2 1 1 1 1 1 2 2 2 1 ...
##  $ stateur: num  4.5 10.5 7.2 5.8 6.5 7.5 5.8 5.8 7.7 6 ...
head(B)
##   sex nwhite age1 age2 tenure1 tenure2 stateur
## 1   1     -1    1    1       1       1     4.5
## 2   1     -1   -1   -1      -1      -1    10.5
## 3  -1     -1    1   -1      -1      -1     7.2
## 4  -1      1    1    1      -1      -1     5.8
## 5   1     -1   -1    1      -1      -1     6.5
## 6   1     -1    1    1       1      -1     7.5

2. Design

Hypothesis and Method

In this experiment, we are to find out the effect of blue collar workers’ sex, color, age, and years of tenure in job lost on the response variable: state unemployment rate (in %). Our hypothesis is that these factors do not affect the response variable. To build an economic experiment, we decide to have only 2 levels for each factors, which results in a \(2^k\) design. To save even more resources, we constructed fractional \(2^k\) design, so we only conduct \(\frac{1}{8}\) of the experimental runs of the full \(2^6\) design, which is a \(2^{6-3}\) design.

Screening

In a fractional factorial experiment, we test only fractions of the full \(2^k\) factorial design. However, there is never free lunch: as we save time and money, we loose the ability to estimate some effects and/or interactions, which is called aliasing. Luckily, not all factors and interactions necessarily have significantly effects our response variable. Thus, in this case where we have to sacrifice some effects, they should be the ones that do not have significant effects. Our goal is to choose the appropriate experimental runs, use the minimum number of experimental runs, and design the experiments so that the effect of factors that have significant effects can be estimated. Therefore, the first step is to identify these “important factors”. This is called a screening process.

##Box plot
par(mfrow=c(1,2))
plot(B$stateur~B$sex, xlab="sex", ylab="unemployment rate")
plot(B$stateur~B$nwhite, xlab="color", ylab="unemployment rate")

par(mfrow=c(1,2))
plot(B$stateur~B$age1, xlab="age1", ylab="unemployment rate")
plot(B$stateur~B$age2, xlab="age2", ylab="unemployment rate")

par(mfrow=c(1,2))
plot(B$stateur~B$tenure1, xlab="tenure1", ylab="unemployment rate")
plot(B$stateur~B$tenure2, xlab="tenure2", ylab="unemployment rate")

##ANOVA
inter = lm(stateur~sex*nwhite+sex*age1+sex*age2+sex*tenure1+sex*tenure2
           +nwhite*age1+nwhite*age2+nwhite*tenure1+nwhite*tenure2
           +age1*age2+age1*tenure1+age1*tenure2
           +age2*tenure1+age2*tenure2
           +tenure1*tenure2, data=B1)
anova(inter)
## Analysis of Variance Table
## 
## Response: stateur
##                  Df  Sum Sq Mean Sq F value   Pr(>F)   
## sex               1    47.9  47.906  7.6778 0.005612 **
## nwhite            1     0.7   0.734  0.1177 0.731597   
## age1              1     1.6   1.597  0.2560 0.612930   
## age2              1    14.5  14.528  2.3284 0.127100   
## tenure1           1    30.1  30.077  4.8203 0.028174 * 
## tenure2           1    15.4  15.447  2.4756 0.115688   
## sex:nwhite        1    20.3  20.340  3.2599 0.071057 . 
## sex:age1          1     0.1   0.084  0.0135 0.907562   
## sex:age2          1     0.6   0.595  0.0954 0.757390   
## sex:tenure1       1     0.1   0.081  0.0130 0.909351   
## sex:tenure2       1     0.1   0.129  0.0207 0.885491   
## nwhite:age1       1     4.0   4.005  0.6419 0.423080   
## nwhite:age2       1     7.0   7.046  1.1292 0.288002   
## nwhite:tenure1    1     4.2   4.225  0.6771 0.410645   
## nwhite:tenure2    1     0.2   0.194  0.0310 0.860157   
## age1:tenure1      1     3.0   2.980  0.4777 0.489519   
## age1:tenure2      1    30.5  30.475  4.8842 0.027150 * 
## age2:tenure1      1     0.0   0.049  0.0078 0.929565   
## age2:tenure2      1     2.9   2.912  0.4667 0.494527   
## Residuals      4857 30305.7   6.240                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the ANOVA results, we notice that main effects of factors: sex, and tenure1are significant. Meanwhile, interaction effect: age1:tenure2 is also significant. Naturally, we want to estimate all these effects, but the truth is, with a \(2^{6-3}\) design, we cannot. To explain this, we need to first understand how fractional factorial designs are generated. A fractional factorial experiment is generated from a full factorial experiment by choosing an alias structure. The alias structure determines which effects are confounded with each other. Since high order effects (>3 factor interactions) are usually considered insignificant and negligible, we sacrifice them in the aliase structure. In a \(2^{6-3}\) design (factor names are A, B, C, D, E, F), the highest possible order interaction effect is ABC, so we first choose \(D=ABC\), \(E=AB\), and \(F=AC\). Clearly, 2-factor interaction effects (AB, AC) are aliased with main effects (E, F), suggesting that the resolution of our design is III. Did chosing a designe generator \(D=ABC\) help increase the resolution? The answer is no. Therefore, to keep everything consistant, we have design generators as \(D=AB\), \(E=AC\), and \("F=BC"\). Defining relation becomes: \(I=ABD=ACE=BCF\). Although the definition of a relolution III fractional factorial design says that main effects are aliased with interaction effects, you may think that we can design the experiment to work around the aliased effects, and only estimate the effects we selected. And that is what we first did: we try to incorporate only the significant effects (main and interaction) in our fractional factorial design, but R retures error suggesting that the number of experimental runs is too small to design an experiment with all these effects estimable.

library(FrF2)
mm <- FrF2(8,6,factor.names = c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2"), 
           estimable=formula("~sex+age1+tenure1+tenure2+age1:tenure2"))

Therefore, we decided to design an experiment focused on the main effects.

mm <- FrF2(8,6,factor.names = c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2"))

We then check the aliasing structure in R, and the result confirms our thoughts on \(2^{6-3}\) design. 2-factor interactions are aliased with each other, and main effects are also aliased with 2 factor interaction effects.

aliasprint(mm)
## $legend
## [1] A=sex     B=nwhite  C=age1    D=age2    E=tenure1 F=tenure2
## 
## $main
## [1] A=BD=CE B=AD=CF C=AE=BF D=AB=EF E=AC=DF F=BC=DE
## 
## $fi2
## [1] AF=BE=CD

Randomize Replicates and Block

To achieve randomization, we randomize the execution order.

mm <- mm[sample(nrow(mm)),]

With the random excution order, we collect data. We first merge the the data (B) with our design (mm). The output is a dataframe containning all data that matches the levels in our design. However, the purpose of fractional factorial deisgn is to save experimental runs. Therefore, we remove all duplicate entrees for the selected conditions, leaving us with 1 experimental run for each condition, which means our experiments have no replicates. Blocking is not used in this experiment.

mmm <- merge(mm, B, by = c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2"), all = FALSE)
##select 1 response variable data for each combination of factors.
unique <- unique( mmm[ , 1:6], nmax = 2)
stateur = mmm$stateur[as.numeric(rownames(unique))]
finalB = cbind(unique,stateur)
finalB
##     sex nwhite age1 age2 tenure1 tenure2 stateur
## 1    -1     -1   -1    1       1       1     8.1
## 16   -1     -1    1    1      -1      -1    11.8
## 183  -1      1   -1   -1       1      -1    11.7
## 193  -1      1    1   -1      -1       1     8.0
## 204   1     -1   -1   -1      -1       1    11.4
## 418   1     -1    1   -1       1      -1    12.3
## 489   1      1   -1    1      -1      -1     5.5
## 532   1      1    1    1       1       1    14.7

3. Analysis

Testing

Using our new dataframe we create a linear model and perform an analysis of variance on it.

model = lm(stateur~sex+nwhite+age1+age2+tenure1+tenure2, data=finalB)
anova(model)
## Analysis of Variance Table
## 
## Response: stateur
##           Df  Sum Sq Mean Sq F value Pr(>F)
## sex        1  2.3112  2.3112  0.0750 0.8298
## nwhite     1  1.7113  1.7113  0.0555 0.8527
## age1       1 12.7512 12.7512  0.4139 0.6361
## age2       1  1.3612  1.3612  0.0442 0.8681
## tenure1    1 12.7513 12.7513  0.4139 0.6361
## tenure2    1  0.1013  0.1013  0.0033 0.9635
## Residuals  1 30.8112 30.8112

As we can see from this new anova we get different results than our original. Here, p-values suggest that sex and tenure1have no significant effects on unemployment rate of blue collar workers. We fail to reject the null hypothesis about factor sex and tenure1 along with the null hypothesis for other factors.

Diagnostics / Model Adequacy Checking

Next we graph Q-Q plots to check our data in the model for normality. If the data is not normal the results of the analysis may not be valid due to invalidation of assumptions made during analysis of variance. We can see below that the Q-Q plot of original model deviates from the linear line, while the fractional factorial design plot is not linear at all. Next we plot a fitted data against the residuals. We can observe a large amount of variation among the plot of original model (residuals at the ends lare narrower than those in the middle). Meanwhile, the plot of the fractional factorial design model look a lot more evenly distributed. Since these results cannot suggest that our fractional factorial design fits the linear model better than the original data, data are more normally distributed in the fractional factorial design. Therefore, ANOVA results of the fractional factorial design model is not more trustworthy.

### for original model
qqnorm(residuals(inter))
qqline(residuals(inter))

### for fractional factorial design
qqnorm(residuals(model))
qqline(residuals(model))

### for original model
plot(fitted(inter),residuals(inter))

### for fractional factorial design
plot(fitted(model),residuals(model))

4. References to Literature

  1. Confounding (also called aliasing). (n.d.). Retrieved December 01, 2016, from http://www.itl.nist.gov/div898/handbook/pri/section3/pri3343.htm

5. Contingencies

Based on the model adequacy testing, it appears that the data is not normally distributed, and that this may not be the best fit of the data. A Kruskal-Wallis test can be used to evaluate whether groups’ distributions are identical, without the assumption that the data is normally distributed. With the kruskal-wallis test, we could only analyze the main effect of one factor at a time. We chose to analyze sex and tenure1.

kruskal.test(stateur~sex, data=finalB)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  stateur by sex
## Kruskal-Wallis chi-squared = 0.33333, df = 1, p-value = 0.5637
kruskal.test(stateur~tenure1, data=finalB)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  stateur by tenure1
## Kruskal-Wallis chi-squared = 2.0833, df = 1, p-value = 0.1489

The p-value for the factor, sex, of this Kruskal-Wallis test is larger than an alpha of 0.05. Therefore, we fail to reject the null, that the variation in blue collar worker unemployment rate is due to error and randomization.

6. Appendices

library(Ecdat)
B <- Ecdat::Benefits

colnames(B)

B <- Ecdat::Benefits [c("sex", "nwhite", "age", "tenure", "stateur")]
head(B) 
str(B)

##transform age into 3-level factor
for (i in 1:4877) {
  if (B$age[i] <= 30) {
  B$age[i] <- "<=30"
  } else if (B$age[i] > 40) {
    B$age[i] <- ">40"
    } else {
       B$age[i] <- ">30 & <=40"
    }
}

##transform tenure into 3-level factor
for (i in 1:4877) {
  if (B$tenure[i] <= 2.5) {
  B$tenure[i] <- "<=2.5"
  } else if (B$tenure[i] > 6) {
    B$tenure[i] <- ">6"
    } else {
       B$tenure[i] <- ">2.5 & <=6"
    }
  }

B1 <- B[c("sex", "nwhite", "age", "tenure","stateur")]
B1$sex <- B$sex
B1$nwhite <- B$nwhite
B1$age <- factor(B$age)
B1$tenure <- factor(B$tenure)
str(B1)
##sex
levels(B1$sex)
##nwhite
levels(B1$nwhite)
##age
levels(B1$age)
##tenure
levels(B1$tenure)

##Transform 3-level factor (age) into 2 2-level factors (age1 and age2)
for (i in 1:4877) {
  if (B$age[i]=="<=30"){
    B$age1[i]<-"30"
    B$age2[i]<-"30"
    } else if (B$age[i]==">40") {
      B$age2[i]<-"40"
      B$age1[i]<-"40"
      } else {
        a <- c(1,2)
        b <- c("30", "40")
        c<- sample(a,1)
        B$age1[i]<- b[c]
        B$age2[i]<- b[3-c]
        }
}

for (i in 1:4877) {
  if (B$tenure[i]=="<=2.5"){
    B$tenure1[i]<-"2.5"
    B$tenure2[i]<-"2.5"
  } else if (B$tenure[i]==">6") {
    B$tenure2[i]<-"6"
    B$tenure1[i]<-"6"
  } else {
    a <- c(1,2)
    b <- c("2.5", "6")
    c<- sample(a,1)
    B$tenure1[i]<- b[c]
    B$tenure2[i]<- b[3-c]
  }
}
B2 <- B[c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2","stateur")]
B2$sex <- B$sex
B2$nwhite <- B$nwhite
B2$age1 <- factor(B$age1)
B2$age2 <- factor(B$age2)
B2$tenure1 <- factor(B$tenure1)
B2$tenure2 <- factor(B$tenure2)
str(B2)

##code levels
B$sex <- as.character(B$sex)
for (i in 1:4877){
  if (as.character(B$sex[i]) == "female") {
    B$sex[i] <- "-1"
  } else {
    B$sex[i] <- "1"
  }
}

B$nwhite <- as.character(B$nwhite)
for (i in 1:4877){
  if (as.character(B$nwhite[i]) == "no") {
    B$nwhite[i] <- "-1"
  } else {
    B$nwhite[i] <- "1"
  }
}

for (i in 1:4877) {
  if (B$age1[i] == "30") {
    B$age1[i] <- "-1"
  } else {
    B$age1[i] <- "1"
  }
}

for (i in 1:4877) {
  if (B$age2[i] == "30") {
    B$age2[i] <- "-1"
  } else {
    B$age2[i] <- "1"
  }
}

for (i in 1:4877) {
  if (B$tenure1[i] =="2.5") {
    B$tenure1[i] <- "-1"
  } else {
    B$tenure1[i] <- "1"
  }
}

for (i in 1:4877) {
  if (B$tenure2[i] =="2.5") {
    B$tenure2[i] <- "-1"
  } else {
    B$tenure2[i] <- "1"
  }
}
B <- B[c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2","stateur")]
B$sex <- factor(B$sex)
B$nwhite <- factor(B$nwhite)
B$age1 <- factor(B$age1)
B$age2 <- factor(B$age2)
B$tenure1 <- factor(B$tenure1)
B$tenure2 <- factor(B$tenure2)
str(B)
head(B)

##Screening
##Box plot
par(mfrow=c(1,2))
plot(B$stateur~B$sex, xlab="sex", ylab="unemployment rate")
plot(B$stateur~B$nwhite, xlab="color", ylab="unemployment rate")
par(mfrow=c(1,2))
plot(B$stateur~B$age1, xlab="age1", ylab="unemployment rate")
plot(B$stateur~B$age2, xlab="age2", ylab="unemployment rate")
par(mfrow=c(1,2))
plot(B$stateur~B$tenure1, xlab="tenure1", ylab="unemployment rate")
plot(B$stateur~B$tenure2, xlab="tenure2", ylab="unemployment rate")
##ANOVA
inter = lm(stateur~sex*nwhite+sex*age1+sex*age2+sex*tenure1+sex*tenure2
           +nwhite*age1+nwhite*age2+nwhite*tenure1+nwhite*tenure2
           +age1*age2+age1*tenure1+age1*tenure2
           +age2*tenure1+age2*tenure2
           +tenure1*tenure2, data=B)
anova(inter)

##fractional factorial design
library(FrF2)
mm <- FrF2(8,6,factor.names = c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2"))
aliasprint(mm)

##randomization
mm <- mm[sample(nrow(mm)),]

##collect data
mmm <- merge(mm, B, by = c("sex", "nwhite", "age1", "age2", "tenure1", "tenure2"), all = FALSE)
nrow(mmm)
##select 1 response variable data for each combination of factors.
unique <- unique( mmm[ , 1:6], nmax = 2)
stateur = mmm$stateur[as.numeric(rownames(unique))]
finalB = cbind(unique,stateur)
finalB

##anova
model = lm(stateur~sex+nwhite+age1+age2+tenure1+tenure2, data=mmm)
anova(model)

### for original model
qqnorm(residuals(inter))
qqline(residuals(inter))
### for fractional factorial design
qqnorm(residuals(model))
qqline(residuals(model))

### for original model
plot(fitted(inter),residuals(inter))
### for fractional factorial design
plot(fitted(model),residuals(model))