# help with unbalanced 2x3x4 factorial ANCOVA

I am attempting an unbalanced 2x3x4 factorial ANCOVA for my toxicology research (contaminant concentrations - y). I have a continuous response variable (blood Hg concentrations - ppm), one continuous covariate (BCI - body condition index), and three categorical independent variables (Sex [male/female], State [OK, NM, TX], and Month [May, June, July, August]). My categorical variables have a different number of levels (Sex = 2, State = 3, and Month = 4). Moreover, I have some empty cell combinations in my analysis from differential sample sizes (ie: no samples in August from one state). I have 97 concentrations in total.

I'm self-taught in R programming, but I do have some codes included below. While I have an initial code framework, I need some help refining my code for interpretation. I can't figure out how to proceed with my analysis in R. I would greatly appreciate any insights you may have!

> datum=read.csv(file.choose())
Blood_Hg State_text State Sex_text Sex Month_text Month BCI
1  0.401   Oklahoma     1   Female   1        May     5   0.3226562
2  0.377   Oklahoma     1   Female   1        May     5   0.3308271
3  0.260   Oklahoma     1   Female   1        May     5   0.3587786
4  0.135   Oklahoma     1     Male   2        May     5   0.3435114
5  0.181   Oklahoma     1   Female   1        May     5   0.3770492
6  0.220   Oklahoma     1     Male   2        May     5   0.3525180
> #Fit maximal model with all parameters for every combination.
> #Blood_Hg as response variable, BCI as continuous covariate, and State, Month, and Sex as categorical independent variables.
> lm1=lm(Blood_Hg~BCI*factor(State)*factor(Month)*factor(Sex), data=datum)
> summary(lm1)

Call:
lm(formula = Blood_Hg ~ BCI * factor(State) * factor(Month) *
factor(Sex), data = datum)
Residuals:
Min       1Q   Median       3Q      Max
-0.19800 -0.03669  0.00000  0.04657  0.48844

Coefficients: (18 not defined because of singularities)
Estimate Std. Error t value  Pr(>|t|)
(Intercept)                                     0.878229   0.369601   2.376  0.0204 *
BCI                                            -1.654220   1.031819  -1.603  0.1136
factor(State)2                                  0.532821   0.676794   0.787  0.4339
factor(State)3                                 -0.291557   2.241489  -0.130  0.8969
factor(Month)6                                 -0.626709   0.469278  -1.335  0.1862
factor(Month)7                                 -0.835299   0.685812  -1.218  0.2275
factor(Month)8                                 -0.800204   0.745027  -1.074  0.2866
factor(Sex)2                                   -0.359534   0.662492  -0.543  0.5891
BCI:factor(State)2                             -2.167151   1.983046  -1.093  0.2784
BCI:factor(State)3                              0.119230   7.143772   0.017  0.9867
BCI:factor(Month)6                              1.718130   1.342873   1.279  0.2052
BCI:factor(Month)7                              2.627388   1.947460   1.349  0.1818
BCI:factor(Month)8                              2.460383   2.101080   1.171  0.2457
factor(State)2:factor(Month)6                  -1.796202   3.318401  -0.541  0.5901
factor(State)3:factor(Month)6                  -0.048460   2.491805  -0.019  0.9845
factor(State)2:factor(Month)7                         NA         NA      NA      NA
factor(State)3:factor(Month)7                   0.457322   2.464728   0.186  0.8534
factor(State)2:factor(Month)8                         NA         NA      NA      NA
factor(State)3:factor(Month)8                         NA         NA      NA      NA
BCI:factor(Sex)2                                0.999650   1.865979   0.536  0.5939
factor(State)2:factor(Sex)2                    -3.305659   1.775924  -1.861  0.0671 .
factor(State)3:factor(Sex)2                     0.439969   1.706550   0.258  0.7973
factor(Month)6:factor(Sex)2                     1.190924   1.158250   1.028  0.3075
factor(Month)7:factor(Sex)2                     0.710682   1.003544   0.708  0.4813
factor(Month)8:factor(Sex)2                    -0.001622   0.123547  -0.013  0.9896
BCI:factor(State)2:factor(Month)6               5.853609  10.426603   0.561  0.5764
BCI:factor(State)3:factor(Month)6               0.177465   7.862417   0.023  0.9821
BCI:factor(State)2:factor(Month)7                     NA         NA      NA      NA
BCI:factor(State)3:factor(Month)7              -1.342813   7.681466  -0.175  0.8618
BCI:factor(State)2:factor(Month)8                     NA         NA      NA      NA
BCI:factor(State)3:factor(Month)8                     NA         NA      NA      NA
BCI:factor(State)2:factor(Sex)2                10.848857   5.777848   1.878  0.0648 .
BCI:factor(State)3:factor(Sex)2                -1.062684   5.144912  -0.207  0.8370
BCI:factor(Month)6:factor(Sex)2                -3.275645   3.296055  -0.994  0.3239
BCI:factor(Month)7:factor(Sex)2                -2.143932   2.824585  -0.759  0.4505
BCI:factor(Month)8:factor(Sex)2                       NA         NA      NA      NA
factor(State)2:factor(Month)6:factor(Sex)2            NA         NA      NA      NA
factor(State)3:factor(Month)6:factor(Sex)2     -0.238484   0.366691  -0.650  0.5177
factor(State)2:factor(Month)7:factor(Sex)2            NA         NA      NA      NA
factor(State)3:factor(Month)7:factor(Sex)2            NA         NA      NA      NA
factor(State)2:factor(Month)8:factor(Sex)2            NA         NA      NA      NA
factor(State)3:factor(Month)8:factor(Sex)2            NA         NA      NA      NA
BCI:factor(State)2:factor(Month)6:factor(Sex)2        NA         NA      NA      NA
BCI:factor(State)3:factor(Month)6:factor(Sex)2        NA         NA      NA      NA
BCI:factor(State)2:factor(Month)7:factor(Sex)2        NA         NA      NA      NA
BCI:factor(State)3:factor(Month)7:factor(Sex)2        NA         NA      NA      NA
BCI:factor(State)2:factor(Month)8:factor(Sex)2        NA         NA      NA      NA
BCI:factor(State)3:factor(Month)8:factor(Sex)2        NA         NA      NA      NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1044 on 67 degrees of freedom
Multiple R-squared:  0.4871,    Adjusted R-squared:  0.2651
F-statistic: 2.194 on 29 and 67 DF,  p-value: 0.004252

> #Note singularities (ie: empty cells).
> #Use step to attempt model simplification.

> lm2=step(lm1)
Start:  AIC=-414.22
Blood_Hg ~ BCI * factor(State) * factor(Month) * factor(Sex)

Step:  AIC=-414.22
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + BCI:factor(Month) + factor(State):factor(Month) +
BCI:factor(Sex) + factor(State):factor(Sex) + factor(Month):factor(Sex) +
BCI:factor(State):factor(Month) + BCI:factor(State):factor(Sex) +
BCI:factor(Month):factor(Sex) + factor(State):factor(Month):factor(Sex)

Df Sum of Sq     RSS     AIC
- BCI:factor(State):factor(Month)          2  0.004065 0.73448 -417.68
- BCI:factor(Month):factor(Sex)            2  0.012669 0.74309 -416.55
- factor(State):factor(Month):factor(Sex)  1  0.004611 0.73503 -415.61
<none>                                                 0.73042 -414.22
- BCI:factor(State):factor(Sex)            2  0.041120 0.77154 -412.91

Step:  AIC=-417.68
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + BCI:factor(Month) + factor(State):factor(Month) +
BCI:factor(Sex) + factor(State):factor(Sex) + factor(Month):factor(Sex) +
BCI:factor(State):factor(Sex) + BCI:factor(Month):factor(Sex) +
factor(State):factor(Month):factor(Sex)

Df Sum of Sq     RSS     AIC
- BCI:factor(Month):factor(Sex)            2  0.012287 0.74677 -420.07
- factor(State):factor(Month):factor(Sex)  1  0.012421 0.74690 -418.05
<none>                                                 0.73448 -417.68
- BCI:factor(State):factor(Sex)            2  0.043754 0.77824 -416.07

Step:  AIC=-420.07
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + BCI:factor(Month) + factor(State):factor(Month) +
BCI:factor(Sex) + factor(State):factor(Sex) + factor(Month):factor(Sex) +
BCI:factor(State):factor(Sex) + factor(State):factor(Month):factor(Sex)

Df Sum of Sq     RSS     AIC
- BCI:factor(Month)                        3  0.022928 0.76970 -423.14
- factor(State):factor(Month):factor(Sex)  2  0.010095 0.75687 -422.77
<none>                                                 0.74677 -420.07
- BCI:factor(State):factor(Sex)            2  0.043285 0.79005 -418.61

Step:  AIC=-423.14
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex)
BCI:factor(State) + factor(State):factor(Month) + BCI:factor(Sex) +
factor(State):factor(Sex) + factor(Month):factor(Sex) + BCI:factor(State):factor(Sex) +
factor(State):factor(Month):factor(Sex)

Df Sum of Sq     RSS     AIC
- factor(State):factor(Month):factor(Sex)  2  0.009007 0.77870 -426.01
<none>                                                 0.76970 -423.14
- BCI:factor(State):factor(Sex)            2  0.042485 0.81218 -421.93

Step:  AIC=-426.01
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + factor(State):factor(Month) + BCI:factor(Sex) +
factor(State):factor(Sex) + factor(Month):factor(Sex) + BCI:factor(State):factor(Sex)

Df Sum of Sq     RSS     AIC
- factor(Month):factor(Sex)      3  0.003481 0.78219 -431.58
- factor(State):factor(Month)    3  0.016751 0.79546 -429.94
<none>                                       0.77870 -426.01
- BCI:factor(State):factor(Sex)  2  0.039874 0.81858 -425.17

Step:  AIC=-431.58
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + factor(State):factor(Month) + BCI:factor(Sex) +
factor(State):factor(Sex) + BCI:factor(State):factor(Sex)

Df Sum of Sq     RSS     AIC
- factor(State):factor(Month)    3  0.014304 0.79649 -435.82
<none>                                       0.78219 -431.58
- BCI:factor(State):factor(Sex)  2  0.041553 0.82374 -430.56

Step:  AIC=-435.82
Blood_Hg ~ BCI + factor(State) + factor(Month) + factor(Sex) +
BCI:factor(State) + BCI:factor(Sex) + factor(State):factor(Sex) +
BCI:factor(State):factor(Sex)

Df Sum of Sq     RSS     AIC
<none>                                       0.79649 -435.82
- BCI:factor(State):factor(Sex)  2  0.041358 0.83785 -434.91
- factor(Month)

> summary(lm2)

Call:
lm(formula = Blood_Hg ~ BCI + factor(State) + factor(Month) +
factor(Sex) + BCI:factor(State) + BCI:factor(Sex) + factor(State):factor(Sex) +
BCI:factor(State):factor(Sex), data = datum)

Residuals:
Min       1Q   Median       3Q      Max
-0.22170 -0.05286  0.00262  0.05873  0.46378

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                      0.359376   0.191356   1.878   0.0639 .
BCI                             -0.209451   0.536244  -0.391   0.6971
factor(State)2                   0.157258   0.377417   0.417   0.6780
factor(State)3                  -0.485631   0.718594  -0.676   0.5011
factor(Month)6                  -0.006954   0.027771  -0.250   0.8029
factor(Month)7                   0.075320   0.034117   2.208   0.0301 *
factor(Month)8                   0.074246   0.038634   1.922   0.0581 .
factor(Sex)2                     0.160179   0.364288   0.440   0.6613
BCI:factor(State)2              -0.956340   1.160217  -0.824   0.4122
BCI:factor(State)3               0.731944   2.212252   0.331   0.7416
BCI:factor(Sex)2                -0.436957   1.028592  -0.425   0.6721
factor(State)2:factor(Sex)2     -3.130893   1.560125  -2.007   0.0481 *
factor(State)3:factor(Sex)2      0.512243   0.932367   0.549   0.5842
BCI:factor(State)2:factor(Sex)2 10.113341   5.175153   1.954   0.0541 .
BCI:factor(State)3:factor(Sex)2 -1.467270   2.836658  -0.517   0.6064
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09856 on 82 degrees of freedom
Multiple R-squared:  0.4407,    Adjusted R-squared:  0.3452
F-statistic: 4.616 on 14 and 82 DF,  p-value: 4.284e-06


Thanks!

• Do you need help just interpreting your results, or with the methodology (since you mentioned unbalanced factorial)? It's not clear what the question here is..
– Jon
Commented Oct 20, 2016 at 0:34
• I could use some help with both. I'm not sure how to interpret the model simplification results and how to then proceed to final results. I'm not familiar with which codes or packages to use in R to move forward. I haven't found any tutorials applicable to my data online. Any help you could provide would be great! Thanks!
– Lou
Commented Oct 20, 2016 at 15:45
• That seems like an awfully big question. Is this your first experience with Design of Experiments, factorial designs?
– Jon
Commented Oct 21, 2016 at 2:48

This may be too big of a question to solve in a single post. You are needing help with

1. Methodology and interpretation
2. Coding

Instead of giving you a long winded response, I'll share some useful references. During my recent MS Statistics program, I skipped Design of Experiments because I disliked the professor (and subject), so I had to self-teach some of the material over the past months.

Here ya go, take the time to actually go through these and not just finish the project as fast as you can. It all has R code which you can use.

Factorial Designs This is a light intro http://www.r-tutor.com/elementary-statistics/analysis-variance/factorial-design

Balanced Factorial Designs https://ww2.coastal.edu/kingw/statistics/R-tutorials/factorial.html

Unbalanced Factorial Designs https://ww2.coastal.edu/kingw/statistics/R-tutorials/unbalanced.html

ANOVA and ANCOVA example http://www.stat.columbia.edu/~martin/W2024/R8.pdf

Now, your model: You're using ANCOVA which is really ANOVA/linear regression with categorical and continuous variables. If it helps for you, dumb it down to regular linear regression to address the relationships between the variables.

I see you included interactions between the covariates. Is this important for your analysis? If so, you should include interaction plots!

Your p-values for your first model are non-significant. Your then run a step-wise regression; if you are not super confident of the mechanics of step-wise regression, I would avoid using it as it can get you in trouble sometimes. If you're comfortable and confident in that last model, lm2, then you can go ahead and start making your conclusions. Otherwise, take the time to actually go through the analysis instead of hitting the fast-forward button.

Anyways, you're missing in your code anova(lm1, lm2) (this will give you a comparison) and anova(lm2) (this will give you the ANOVA table). If you took the time to go through the references I gave you, you shouldn't have a problem figuring out what it all means.

Good luck. Effort goes a long way...