0
$\begingroup$

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())
> head(datum)
  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!

$\endgroup$
3
  • $\begingroup$ 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.. $\endgroup$
    – Jon
    Commented Oct 20, 2016 at 0:34
  • $\begingroup$ 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! $\endgroup$
    – Lou
    Commented Oct 20, 2016 at 15:45
  • $\begingroup$ That seems like an awfully big question. Is this your first experience with Design of Experiments, factorial designs? $\endgroup$
    – Jon
    Commented Oct 21, 2016 at 2:48

1 Answer 1

2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.