HW 10

A| Millipedes

• Test the research hypothesis that species affects alanine concentrations.

• Test the research hypothesis that sex affects alanine concentrations.

• Test the research hypothesis that the effect of sex depends on species.

A.1| Data Read

millipedes <- read.csv("Data/millipede.csv")
str(millipedes) #great!

## 'data.frame':    24 obs. of  4 variables:
##  $ EU       : chr  "A" "B" "C" "D" ...
##  $ species  : chr  "alpha" "alpha" "alpha" "alpha" ...
##  $ hemolymph: num  21.5 19.6 20.9 22.8 14.8 15.6 13.5 16.4 14.5 17.4 ...
##  $ sex      : chr  "male" "male" "male" "male" ...

A.2| Assumptions Check

#fit model including species, sex, and their interaction "depends on"
millipedeModel <- lm(hemolymph~species*sex, data=millipedes)

#check assumptions
par(mfrow=c(2,2))
plot(millipedeModel)#tails look a little funky, and there is one or two groups with some outlier residuals

#double check 
shapiro.test(millipedeModel$residuals) #normal yes!

## 
##  Shapiro-Wilk normality test
## 
## data:  millipedeModel$residuals
## W = 0.94757, p-value = 0.2399

bartlett.test(hemolymph ~ interaction(species, sex), data=millipedes) #no significant evidence of unequal variances

## 
##  Bartlett test of homogeneity of variances
## 
## data:  hemolymph by interaction(species, sex)
## Bartlett's K-squared = 0.81732, df = 5, p-value = 0.9759

Assumptions for proceeding with the model have all been met.

A.3| LM

summary(millipedeModel) #tells R2 effect and shows estimates

## 
## Call:
## lm(formula = hemolymph ~ species * sex, data = millipedes)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.525 -1.200  0.050  1.038  1.825 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          15.0750     0.7267  20.746 5.11e-14 ***
## speciesbeta          -2.4000     1.0276  -2.335   0.0313 *  
## specieszeta          -1.3500     1.0276  -1.314   0.2054    
## sexmale               6.1250     1.0276   5.960 1.22e-05 ***
## speciesbeta:sexmale  -2.6250     1.4533  -1.806   0.0876 .  
## specieszeta:sexmale  -1.3250     1.4533  -0.912   0.3740    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.453 on 18 degrees of freedom
## Multiple R-squared:  0.8409, Adjusted R-squared:  0.7967 
## F-statistic: 19.02 on 5 and 18 DF,  p-value: 1.281e-06

#look at Anova to answer research qs
Anova(millipedeModel)#this shows that the interaction is not significant so we can take it out of teh model. the df for species is 2,18 and df for sex is 1,18

## Anova Table (Type II tests)
## 
## Response: hemolymph
##              Sum Sq Df F value    Pr(>F)    
## species      55.261  2 13.0821 0.0003103 ***
## sex         138.720  1 65.6794 2.037e-07 ***
## species:sex   6.891  2  1.6313 0.2233107    
## Residuals    38.018 18                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

millipedeModel2 <- aov(hemolymph ~ species + sex, data = millipedes) #here we use aov to prepare for tukey and take out the interaction term since Anova showed it wasnt significant
summary(millipedeModel2)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## species      2  55.26   27.63   12.30 0.000328 ***
## sex          1 138.72  138.72   61.78 1.53e-07 ***
## Residuals   20  44.91    2.25                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(millipedeModel2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = hemolymph ~ species + sex, data = millipedes)
## 
## $species
##               diff        lwr        upr     p adj
## beta-alpha -3.7125 -5.6080527 -1.8169473 0.0002160
## zeta-alpha -2.0125 -3.9080527 -0.1169473 0.0362524
## zeta-beta   1.7000 -0.1955527  3.5955527 0.0838160
## 
## $sex
##                 diff      lwr      upr p adj
## male-female 4.808333 3.532249 6.084418 2e-07

Species effect: Significant effect by species on alanine concentration, \(F_{2,18}\) = 13.08, p < .001. At least one species mean differs. After running Tukey posthoc, it showed that beta-alpha and zeta-alpha were statistically significant and had different enough means.

Sex effect: Significant effect by sex on alanine concentration, \(F_{1,18}\) = 65.68, p < .001.

Does the effect of sex depend on species? No, the interaction term was not significant after running Anova(millipedeModel) it showed that there was not a reason to include the interaction in the final model.

A.4| Figures

#make letters so readability is best for stats
speciesletters <- data.frame(species = c("alpha", "beta", "zeta"),label = c("a", "b", "b"), height=c(23.1, 18.1, 20.6))
sexletters <- data.frame(sex = c("female", "male"),label = c("b", "a"), height=c(16.8, 23.1))

yrange <- range(millipedes$hemolymph) #share y axis

#species
speciesplot <- ggplot(millipedes, aes(species, hemolymph, fill = species)) +
  geom_boxplot() +
  labs(x = "Species", y = "Alanine concentration (mg/100 ml)") +
  theme_classic() +
  theme(legend.position = "none")+
  coord_cartesian(ylim = yrange) + #get em to share y axis
  scale_fill_manual(values=c("#97c8b3", "#9a890e", "#c59ec0"))+
  geom_text(data = speciesletters,aes(x = species, y = height, label = label)) #add on letters above plot

#sex
sexplot <- ggplot(millipedes, aes(sex, hemolymph, fill = sex)) +
  geom_boxplot() +
  labs(x = "Sex", y = "Alanine concentration (mg/100 ml)") +
  theme_classic() +
  theme(legend.position = "none")+
  coord_cartesian(ylim = yrange)+
  scale_fill_manual(values=c("#e98891","#4368be"))+
  geom_text(data = sexletters,aes(x = sex,y=height, label = label))

(speciesplot|sexplot)+
  plot_layout(axes = "collect")

Figure 1. Alanine concentrations (mg/100 ml) by millipede species and sex. Alanine concentration differed significantly among species, \(F_{2,18}\) = 13.08, p < .001, and between sexes, \(F_{1,18}\) = 65.68, p < .001, with males having higher concentrations than females. The species × sex interaction was not significant, \(F_{2,18}\) = 1.63, p = .223.

B| Blocking & Covariates

Blocks and covariates are another way to reduce \(SS_{error}\) by accounting for variation created by categorical IV (blocks) or meristic IV (covariates). In general factors that are blocks or covariates are not of central interest in your study. Block and covariates ARE NOT included in interaction terms.

These are data from an experiment on hairy mole rats, which are solitary fossorial creatures. Six mole rats were given testosterone implants, six were given control implants. The field metabolic rate of each of the 12 mole rats was measured over a 3 day period. You made two different experimental runs. Because they live in a number of habitat types, the soil temperatures experienced varies from rat to rat. You, of course, took pains to take soil temperatures. Determine if there was an effect of hormone treatment on field metabolic rate (MR). Also see if you should include Average temperature and/or experimental run as blocks or covariates in your model.

B.1| Data

Treatment <- factor(rep(c("Testosterone", "Control"), each = 6))

FieldMR <- c(11.1,11.7,12.5,12.9,13.1,11.2,12.2,10.7,11.4,11,11.2,9.9)

Temp <- c(19.1,22.7,23.8,22.6,24.8,25.1,21.2,21.5,23.9,25.5,26,28.7)

ExperimentNumber <- factor(rep(c(1, 2), times = 6))

molerat <- data.frame(Treatment, FieldMR, Temp, ExperimentNumber)

B.2| Assumptions Check

#fit model as fixed effects and check assumption s
treatmentonlyModel <- lm(FieldMR~Treatment, data=molerat)
treatmentTempModel <- lm(FieldMR~Treatment+Temp, data=molerat)
treatmentIDModel <- lm(FieldMR~Treatment+ExperimentNumber, data=molerat)

par(mfrow = c(2, 2))
plot(treatmentonlyModel)

plot(treatmentTempModel)

plot(treatmentIDModel)

Assumptions have been met and I will proceed with comparing models.

B.3| Model Comparison

Comparing F, p, and MSE to determine best model.

Anova(treatmentonlyModel)

## Anova Table (Type II tests)
## 
## Response: FieldMR
##           Sum Sq Df F value  Pr(>F)  
## Treatment 3.1008  1  4.6408 0.05666 .
## Residuals 6.6817 10                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(treatmentTempModel)

## Anova Table (Type II tests)
## 
## Response: FieldMR
##           Sum Sq Df F value Pr(>F)
## Treatment 2.3879  1  3.3199 0.1018
## Temp      0.2082  1  0.2894 0.6037
## Residuals 6.4735  9

Anova(treatmentIDModel)

## Anova Table (Type II tests)
## 
## Response: FieldMR
##                  Sum Sq Df F value  Pr(>F)  
## Treatment        3.1008  1  5.2847 0.04709 *
## ExperimentNumber 1.4008  1  2.3874 0.15672  
## Residuals        5.2808  9                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sigma(treatmentonlyModel)^2

## [1] 0.6681667

sigma(treatmentTempModel)^2

## [1] 0.7192788

sigma(treatmentIDModel)^2

## [1] 0.5867593

I selected the model including treatment and experimental run. Experimental number accounts for the experimental blocking structure, and including it reduced MSE from 0.6682 to 0.5868. Adding temperature instead increased MSE to 0.7193, so I did not retain temperature. After accounting for experimental number, testosterone-treated rats had a mean field metabolic rate approximately 1.02 units higher than controls, F(1, 9) = 5.28, p = .047. Therefore, I rejected the null hypothesis of equal treatment means at α = .05.

Final Model

treatmentIDModel <- lm(FieldMR~Treatment+ExperimentNumber, data=molerat)

C| ANCOVA

Run a linear model with an interaction between a categorical and numerical predictor (ANCOVA) on the effects of water depth and whether a not the stream is in a burned area. I put the code and the data are below. You can also look at the worksheet attached to the ANCOVA lecture.

fish<-read.csv("Data/ancovafish.csv")

str(fish)

## 'data.frame':    22 obs. of  3 variables:
##  $ water: int  0 1 2 3 4 5 6 7 8 9 ...
##  $ burn : chr  "unburn" "unburn" "unburn" "unburn" ...
##  $ Fish : num  37 35.9 38.2 40 41.1 46.5 48 50.1 51.1 55.2 ...

modfish<-lm(Fish~burn*water, data=fish) #interaction

par(mfrow=c(2,2))
plot(modfish)

summary(modfish)

## 
## Call:
## lm(formula = Fish ~ burn * water, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1091 -0.8045 -0.2173  0.9614  2.4455 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       43.7136     0.7230  60.462  < 2e-16 ***
## burnunburn        -9.1591     1.0225  -8.958 4.72e-08 ***
## water              1.1900     0.1222   9.738 1.34e-08 ***
## burnunburn:water   0.9736     0.1728   5.634 2.41e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.282 on 18 degrees of freedom
## Multiple R-squared:  0.9631, Adjusted R-squared:  0.957 
## F-statistic: 156.6 on 3 and 18 DF,  p-value: 4.386e-13

Anova(modfish, type="III")

## Anova Table (Type III tests)
## 
## Response: Fish
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 6005.6  1 3655.676 < 2.2e-16 ***
## burn         131.8  1   80.243 4.723e-08 ***
## water        155.8  1   94.819 1.343e-08 ***
## burn:water    52.1  1   31.737 2.407e-05 ***
## Residuals     29.6 18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#significant interaction so run the models separately to get parameter estimates

bfish<-subset(fish, burn =="burn")

ubfish<-subset(fish, burn =="unburn")

modburnfish<-lm(Fish~water, data=bfish)

summary(modburnfish)

## 
## Call:
## lm(formula = Fish ~ water, data = bfish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5136 -0.7986 -0.2436  0.7414  1.8464 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  43.7136     0.6851   63.81 2.88e-13 ***
## water         1.1900     0.1158   10.28 2.85e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.214 on 9 degrees of freedom
## Multiple R-squared:  0.9215, Adjusted R-squared:  0.9127 
## F-statistic: 105.6 on 1 and 9 DF,  p-value: 2.85e-06

modunburnfish<-lm(Fish~water, data=ubfish)

summary(modunburnfish)

## 
## Call:
## lm(formula = Fish ~ water, data = ubfish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1091 -0.7909 -0.1909  0.7954  2.4455 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.5545     0.7590   45.52 5.96e-12 ***
## water         2.1636     0.1283   16.86 4.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.346 on 9 degrees of freedom
## Multiple R-squared:  0.9693, Adjusted R-squared:  0.9659 
## F-statistic: 284.4 on 1 and 9 DF,  p-value: 4.063e-08

The relationship between water depth and fish density depended on burn status, F(1, 18) = 31.74, p < .001. Fish density increased significantly with depth in both burned streams, F(1, 9) = 105.6, p < .001, and unburned streams, F(1, 9) = 284.4, p < .001. However, the increase was steeper in unburned streams: predicted fish density increased by 2.16 units per cm of depth, compared with 1.19 units per cm in burned streams. The full model explained 96.3% of the variation in fish density.

Burned: Predicted fish density = 43.71 + 1.19(water depth).

Unburned: Predicted fish density = 34.55 + 2.16(water depth).