OCN 683 - Homework 1
Sarah Pennington
###(1) Start with some data exploration. How many observations are there of each bird species in each experimental treatment? Make a boxplot or violinplot to show the distribution of foraging heights of each bird species in each experimental treatment.What have you learned so far?**
getwd()## [1] "C:/Users/sarah/Desktop/Masters program/Classes_spring_25"bird <- read_csv("foraging.height.edit.csv",
col_types = cols(...1 = col_skip(), Area_ha = col_number(),
Elevation_m = col_number(), tot.arth.bm = col_number(),
foraging.ht.m = col_number()))## New names:
## • `` -> `...1`names(bird) #names of columns## [1] "SPECIES" "Kipuka" "Area_ha" "Elevation_m"
## [5] "Rat_Removal" "DATE" "OBSERVER" "Year"
## [9] "tot.arth.bm" "dietary.grouping" "foraging.ht.m"str(bird) #how many observations, and what are they for each column ## spc_tbl_ [527 × 11] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ SPECIES : chr [1:527] "JAWE" "HAAM" "IIWI" "HAAM" ...
## $ Kipuka : chr [1:527] "K13" "K13" "K5" "K5" ...
## $ Area_ha : num [1:527] 0.1 0.1 10.5 10.5 10.5 ...
## $ Elevation_m : num [1:527] 1575 1575 1632 1632 1632 ...
## $ Rat_Removal : chr [1:527] "treated" "treated" "treated" "treated" ...
## $ DATE : chr [1:527] "1/28/12" "1/23/12" "1/22/12" "1/22/12" ...
## $ OBSERVER : chr [1:527] "JK" "NF" "WK" "WK" ...
## $ Year : num [1:527] 2012 2012 2012 2012 2012 ...
## $ tot.arth.bm : num [1:527] 1.16 1.16 149.63 149.63 149.63 ...
## $ dietary.grouping: chr [1:527] "frugivore" "insectivore" "nectarivore" "insectivore" ...
## $ foraging.ht.m : num [1:527] 2.24 2.74 4.15 4.65 5.15 ...
## - attr(*, "spec")=
## .. cols(
## .. ...1 = col_skip(),
## .. SPECIES = col_character(),
## .. Kipuka = col_character(),
## .. Area_ha = col_number(),
## .. Elevation_m = col_number(),
## .. Rat_Removal = col_character(),
## .. DATE = col_character(),
## .. OBSERVER = col_character(),
## .. Year = col_double(),
## .. tot.arth.bm = col_number(),
## .. dietary.grouping = col_character(),
## .. foraging.ht.m = col_number()
## .. )
## - attr(*, "problems")=<externalptr>dim(bird) #what are the dimensions of the data ## [1] 527 11##How many observations for each bird species in each treatment group?
sum.bird <- bird %>%
group_by(SPECIES, Rat_Removal) %>%
summarise(n = n()) ## `summarise()` has grouped output by 'SPECIES'. You can override using the
## `.groups` argument.print(sum.bird)## # A tibble: 12 × 3
## # Groups: SPECIES [6]
## SPECIES Rat_Removal n
## <chr> <chr> <int>
## 1 APAP treated 131
## 2 APAP untreated 188
## 3 HAAM treated 47
## 4 HAAM untreated 42
## 5 HAEL treated 2
## 6 HAEL untreated 3
## 7 IIWI treated 30
## 8 IIWI untreated 43
## 9 JAWE treated 16
## 10 JAWE untreated 9
## 11 OMAO treated 7
## 12 OMAO untreated 9#Make a boxplot or a violin plot:
plot <- ggplot(bird, aes(x=Rat_Removal,y=foraging.ht.m,fill= SPECIES)) +geom_violin(trim = FALSE ) +theme_clean() + stat_summary(fun=median, geom="point",size=2, position = position_dodge(width = 0.9))+labs(y="Foraging height",
x="treatment",
title="Foraging height by treatment")
print(plot)
#What have we learned so far?
## The foraging height based on bird species, when looking at the spread of the data with a violin plot, doesn't appear to drastically differ between treated and untreated groups. (treated = rats removed). However, we can look at shifts for each individual species when rats are present to see if there are within species trends in foraging height. ###(2) Now make a plot showing the mean foraging height of each bird species in each treatment, and include error bars displaying +/- one standard error of the mean. What is the meaning of the standard error of the mean? How does this plot differ from the plot in #1?
sum.bird <- bird %>%
group_by(SPECIES, Rat_Removal) %>%
summarise(n = n(), mean = mean(foraging.ht.m, na.rm = TRUE),
se = sd(foraging.ht.m, na.rm = TRUE) / sqrt(n())) ## `summarise()` has grouped output by 'SPECIES'. You can override using the
## `.groups` argument.print(sum.bird)## # A tibble: 12 × 5
## # Groups: SPECIES [6]
## SPECIES Rat_Removal n mean se
## <chr> <chr> <int> <dbl> <dbl>
## 1 APAP treated 131 15.7 0.436
## 2 APAP untreated 188 17.3 0.336
## 3 HAAM treated 47 14.8 0.697
## 4 HAAM untreated 42 16.5 0.803
## 5 HAEL treated 2 15.9 2.30
## 6 HAEL untreated 3 14.3 1.34
## 7 IIWI treated 30 17.5 0.777
## 8 IIWI untreated 43 19.3 0.534
## 9 JAWE treated 16 15.6 1.70
## 10 JAWE untreated 9 19.6 1.16
## 11 OMAO treated 7 14.9 1.44
## 12 OMAO untreated 9 17.2 0.954##There are only 5 total observations of species HAEL
plot1 <- ggplot(sum.bird, aes(x = Rat_Removal, y = mean, group = SPECIES, fill= SPECIES )) +
geom_col(position = position_dodge(width = 0.9)) +
geom_errorbar(aes(ymin = mean-se, ymax = mean+se),
position = position_dodge(width = 0.9), width = 0.2)+labs(
y = "Mean Foraging Height (m)",
x = "Treatment",
title = "Mean Foraging Height by Bird Species and Treatment") +
theme_clean()
print(plot1)
#What is the meaning of standard error of the mean?
### Assuming the data are normaly distributed the standard error is the the average amount each datapoint differs from the sample mean (standard deviation) divided by the square root of the number of samples (n). So, when there are more samples, the SE will be smaller.
#How does this plot differ from the first plot?
### The first plot shows the overall spread of the raw data, this plot shows the average foraging height, plus the spread of the data. The information in this plot (means) can be used to figure out if there is a significant difference in foraging height between any of the groups. (3) Fit a linear model that tests whether the bird species differ in foraging height,whether the rat removal treatment affects foraging height, and whether the effect of the rat removal differs between the bird species. Make plots to assess the distribution of the residuals. Report F-tests for the terms of the model. Create an effects plot displaying the fitted effects. How do you interpret the results? What are the magnitudes of the effects?
#fit the linear model
bird.model = lm( foraging.ht.m ~ SPECIES * Rat_Removal, data = bird)
summary(bird.model)##
## Call:
## lm(formula = foraging.ht.m ~ SPECIES * Rat_Removal, data = bird)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.3712 -3.1989 0.4048 3.4124 9.9124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.7426 0.4108 38.326 <2e-16 ***
## SPECIESHAAM -0.9253 0.7994 -1.157 0.2476
## SPECIESHAEL 0.1124 3.3496 0.034 0.9732
## SPECIESIIWI 1.7836 0.9516 1.874 0.0614 .
## SPECIESJAWE -0.1713 1.2450 -0.138 0.8906
## SPECIESOMAO -0.8419 1.8238 -0.462 0.6446
## Rat_Removaluntreated 1.5426 0.5351 2.883 0.0041 **
## SPECIESHAAM:Rat_Removaluntreated 0.1377 1.1326 0.122 0.9033
## SPECIESHAEL:Rat_Removaluntreated -3.0558 4.3249 -0.707 0.4802
## SPECIESIIWI:Rat_Removaluntreated 0.2512 1.2398 0.203 0.8395
## SPECIESJAWE:Rat_Removaluntreated 2.4428 2.0306 1.203 0.2295
## SPECIESOMAO:Rat_Removaluntreated 0.7522 2.4289 0.310 0.7569
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.701 on 515 degrees of freedom
## Multiple R-squared: 0.06603, Adjusted R-squared: 0.04608
## F-statistic: 3.31 on 11 and 515 DF, p-value: 0.0002061#make plots to assess the distribution of the residuals
library(ggResidpanel)
resid_panel(bird.model, plots = c('resid', 'qq', 'lev', 'hist'))
#Report F tests for terms of the model.
Anova(bird.model)## Anova Table (Type II tests)
##
## Response: foraging.ht.m
## Sum Sq Df F value Pr(>F)
## SPECIES 357.3 5 3.2329 0.006952 **
## Rat_Removal 372.2 1 16.8403 4.726e-05 ***
## SPECIES:Rat_Removal 45.9 5 0.4157 0.837908
## Residuals 11382.8 515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##F values
##SPECIES 3.32329 **
##Rat_Removal 16.8403 ***
##SPECIES:Rat_Removal 0.4157
#Create an effects plot displaying fitted effects
library(ggeffects)
plot(ggeffect(bird.model, terms = c('SPECIES', 'Rat_Removal')))
#How do you interpret these results? What are the magnitudes of the effects?
## In general, Rat_removal had an effect of greater magnitude (F = 16.84) than species (F = 3.323) on predicting foraging height. Notably, species HAEL had the largest foraging height predicted values variation. There was only 5 total observations for this species, so that makes sense. (4) So far we have examined the effect of the rat removal, but we have not considered the role of kipuka size, or the role of food availability (which is influenced in part by kipuka size).Make a new model that adds kipuka area and arthopod biomass as predictors, and include interactions that test whether the treatment effect depends on these additional predictors.
###Kipuka area and arthopod biomass are both continuous predictors – before you add them to the model, make some scatterplots to assess whether these predictors should be transformed when including them in the model. (When assessing predictors for transformation, it doesn’t matter if the predictors are normally distributed (this only matters for the response variable) – what matters is whether the predictors are very skewed, such that a few outlying points will have a large influence on a fitted regression line. For skewed predictors, a log or square root transformation will generally help.)
###Report F-tests and effects plots. How do the results of this model differ from the model in #3? How do you interpret the results at this stage?
#Make some scatterplots to see if Kipuku area and arthropod biomass should be transformed.
##kipuku area
kip.plot <- ggplot(bird, aes(x = Area_ha, y = foraging.ht.m)) +geom_point() +
theme_clean()
print(kip.plot)###No transformation needed
## Sarthropod biomass
bio.plot <- ggplot(bird, aes(x = tot.arth.bm, y = foraging.ht.m)) +geom_point() +
theme_clean()
print(bio.plot) ###Should transform
bio.plot <- ggplot(bird, aes(x = log10(tot.arth.bm), y = foraging.ht.m)) +geom_point() +
theme_clean()
print(bio.plot)
#Make a new model including kipuku area and biomass
bird.model.2 <- lm(foraging.ht.m ~ SPECIES * Rat_Removal * Area_ha * log10(tot.arth.bm), data = bird)
summary(bird.model.2)##
## Call:
## lm(formula = foraging.ht.m ~ SPECIES * Rat_Removal * Area_ha *
## log10(tot.arth.bm), data = bird)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.6127 -2.5738 0.0041 2.7795 11.0368
##
## Coefficients: (3 not defined because of singularities)
## Estimate
## (Intercept) 12.72727
## SPECIESHAAM -1.50235
## SPECIESHAEL -103.03837
## SPECIESIIWI 1.00871
## SPECIESJAWE 0.68699
## SPECIESOMAO -1.73937
## Rat_Removaluntreated 1.60422
## Area_ha 0.83686
## log10(tot.arth.bm) 4.88008
## SPECIESHAAM:Rat_Removaluntreated -0.11689
## SPECIESHAEL:Rat_Removaluntreated 115.05665
## SPECIESIIWI:Rat_Removaluntreated 3.68893
## SPECIESJAWE:Rat_Removaluntreated 4.41517
## SPECIESOMAO:Rat_Removaluntreated 1.18937
## SPECIESHAAM:Area_ha 0.13928
## SPECIESHAEL:Area_ha 6.90875
## SPECIESIIWI:Area_ha -0.06974
## SPECIESJAWE:Area_ha -0.96353
## SPECIESOMAO:Area_ha 2.40704
## Rat_Removaluntreated:Area_ha -0.39877
## SPECIESHAAM:log10(tot.arth.bm) 0.54227
## SPECIESHAEL:log10(tot.arth.bm) 94.74617
## SPECIESIIWI:log10(tot.arth.bm) 1.55913
## SPECIESJAWE:log10(tot.arth.bm) -9.99729
## SPECIESOMAO:log10(tot.arth.bm) -14.14989
## Rat_Removaluntreated:log10(tot.arth.bm) -4.80826
## Area_ha:log10(tot.arth.bm) -1.13184
## SPECIESHAAM:Rat_Removaluntreated:Area_ha -0.28832
## SPECIESHAEL:Rat_Removaluntreated:Area_ha -10.43675
## SPECIESIIWI:Rat_Removaluntreated:Area_ha -0.69805
## SPECIESJAWE:Rat_Removaluntreated:Area_ha 0.07669
## SPECIESOMAO:Rat_Removaluntreated:Area_ha -2.23746
## SPECIESHAAM:Rat_Removaluntreated:log10(tot.arth.bm) -2.88334
## SPECIESHAEL:Rat_Removaluntreated:log10(tot.arth.bm) NA
## SPECIESIIWI:Rat_Removaluntreated:log10(tot.arth.bm) -2.46744
## SPECIESJAWE:Rat_Removaluntreated:log10(tot.arth.bm) 11.86150
## SPECIESOMAO:Rat_Removaluntreated:log10(tot.arth.bm) 13.65459
## SPECIESHAAM:Area_ha:log10(tot.arth.bm) -0.11535
## SPECIESHAEL:Area_ha:log10(tot.arth.bm) NA
## SPECIESIIWI:Area_ha:log10(tot.arth.bm) -0.27442
## SPECIESJAWE:Area_ha:log10(tot.arth.bm) 4.72574
## SPECIESOMAO:Area_ha:log10(tot.arth.bm) 3.92278
## Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 1.31802
## SPECIESHAAM:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.47887
## SPECIESHAEL:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) NA
## SPECIESIIWI:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.57888
## SPECIESJAWE:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) -4.43581
## SPECIESOMAO:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) -4.26608
## Std. Error t value
## (Intercept) 0.51028 24.942
## SPECIESHAAM 1.12958 -1.330
## SPECIESHAEL 132.20209 -0.779
## SPECIESIIWI 2.05674 0.490
## SPECIESJAWE 1.74324 0.394
## SPECIESOMAO 3.28008 -0.530
## Rat_Removaluntreated 0.66451 2.414
## Area_ha 0.11191 7.478
## log10(tot.arth.bm) 0.77211 6.320
## SPECIESHAAM:Rat_Removaluntreated 1.59670 -0.073
## SPECIESHAEL:Rat_Removaluntreated 149.86334 0.768
## SPECIESIIWI:Rat_Removaluntreated 2.71620 1.358
## SPECIESJAWE:Rat_Removaluntreated 4.72974 0.933
## SPECIESOMAO:Rat_Removaluntreated 4.23888 0.281
## SPECIESHAAM:Area_ha 0.21514 0.647
## SPECIESHAEL:Area_ha 8.90883 0.775
## SPECIESIIWI:Area_ha 0.29466 -0.237
## SPECIESJAWE:Area_ha 0.47877 -2.012
## SPECIESOMAO:Area_ha 2.22744 1.081
## Rat_Removaluntreated:Area_ha 0.14644 -2.723
## SPECIESHAAM:log10(tot.arth.bm) 1.96729 0.276
## SPECIESHAEL:log10(tot.arth.bm) 125.87141 0.753
## SPECIESIIWI:log10(tot.arth.bm) 2.36067 0.660
## SPECIESJAWE:log10(tot.arth.bm) 2.77309 -3.605
## SPECIESOMAO:log10(tot.arth.bm) 14.10617 -1.003
## Rat_Removaluntreated:log10(tot.arth.bm) 0.89441 -5.376
## Area_ha:log10(tot.arth.bm) 0.19633 -5.765
## SPECIESHAAM:Rat_Removaluntreated:Area_ha 0.27669 -1.042
## SPECIESHAEL:Rat_Removaluntreated:Area_ha 13.30167 -0.785
## SPECIESIIWI:Rat_Removaluntreated:Area_ha 0.37681 -1.853
## SPECIESJAWE:Rat_Removaluntreated:Area_ha 1.16210 0.066
## SPECIESOMAO:Rat_Removaluntreated:Area_ha 2.36914 -0.944
## SPECIESHAAM:Rat_Removaluntreated:log10(tot.arth.bm) 2.49852 -1.154
## SPECIESHAEL:Rat_Removaluntreated:log10(tot.arth.bm) NA NA
## SPECIESIIWI:Rat_Removaluntreated:log10(tot.arth.bm) 3.33171 -0.741
## SPECIESJAWE:Rat_Removaluntreated:log10(tot.arth.bm) 4.51134 2.629
## SPECIESOMAO:Rat_Removaluntreated:log10(tot.arth.bm) 14.26506 0.957
## SPECIESHAAM:Area_ha:log10(tot.arth.bm) 0.34162 -0.338
## SPECIESHAEL:Area_ha:log10(tot.arth.bm) NA NA
## SPECIESIIWI:Area_ha:log10(tot.arth.bm) 0.37579 -0.730
## SPECIESJAWE:Area_ha:log10(tot.arth.bm) 1.24468 3.797
## SPECIESOMAO:Area_ha:log10(tot.arth.bm) 7.41901 0.529
## Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.21657 6.086
## SPECIESHAAM:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.38984 1.228
## SPECIESHAEL:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) NA NA
## SPECIESIIWI:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.46056 1.257
## SPECIESJAWE:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 1.45691 -3.045
## SPECIESOMAO:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 7.46620 -0.571
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## SPECIESHAAM 0.184144
## SPECIESHAEL 0.436126
## SPECIESIIWI 0.624046
## SPECIESJAWE 0.693691
## SPECIESOMAO 0.596160
## Rat_Removaluntreated 0.016145 *
## Area_ha 3.60e-13 ***
## log10(tot.arth.bm) 5.94e-10 ***
## SPECIESHAAM:Rat_Removaluntreated 0.941671
## SPECIESHAEL:Rat_Removaluntreated 0.443015
## SPECIESIIWI:Rat_Removaluntreated 0.175061
## SPECIESJAWE:Rat_Removaluntreated 0.351034
## SPECIESOMAO:Rat_Removaluntreated 0.779149
## SPECIESHAAM:Area_ha 0.517665
## SPECIESHAEL:Area_ha 0.438428
## SPECIESIIWI:Area_ha 0.813018
## SPECIESJAWE:Area_ha 0.044723 *
## SPECIESOMAO:Area_ha 0.280401
## Rat_Removaluntreated:Area_ha 0.006701 **
## SPECIESHAAM:log10(tot.arth.bm) 0.782941
## SPECIESHAEL:log10(tot.arth.bm) 0.451984
## SPECIESIIWI:log10(tot.arth.bm) 0.509274
## SPECIESJAWE:log10(tot.arth.bm) 0.000345 ***
## SPECIESOMAO:log10(tot.arth.bm) 0.316316
## Rat_Removaluntreated:log10(tot.arth.bm) 1.19e-07 ***
## Area_ha:log10(tot.arth.bm) 1.46e-08 ***
## SPECIESHAAM:Rat_Removaluntreated:Area_ha 0.297906
## SPECIESHAEL:Rat_Removaluntreated:Area_ha 0.433062
## SPECIESIIWI:Rat_Removaluntreated:Area_ha 0.064561 .
## SPECIESJAWE:Rat_Removaluntreated:Area_ha 0.947411
## SPECIESOMAO:Rat_Removaluntreated:Area_ha 0.345429
## SPECIESHAAM:Rat_Removaluntreated:log10(tot.arth.bm) 0.249063
## SPECIESHAEL:Rat_Removaluntreated:log10(tot.arth.bm) NA
## SPECIESIIWI:Rat_Removaluntreated:log10(tot.arth.bm) 0.459302
## SPECIESJAWE:Rat_Removaluntreated:log10(tot.arth.bm) 0.008830 **
## SPECIESOMAO:Rat_Removaluntreated:log10(tot.arth.bm) 0.338944
## SPECIESHAAM:Area_ha:log10(tot.arth.bm) 0.735762
## SPECIESHAEL:Area_ha:log10(tot.arth.bm) NA
## SPECIESIIWI:Area_ha:log10(tot.arth.bm) 0.465584
## SPECIESJAWE:Area_ha:log10(tot.arth.bm) 0.000165 ***
## SPECIESOMAO:Area_ha:log10(tot.arth.bm) 0.597224
## Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 2.37e-09 ***
## SPECIESHAAM:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.219907
## SPECIESHAEL:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) NA
## SPECIESIIWI:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.209398
## SPECIESJAWE:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.002457 **
## SPECIESOMAO:Rat_Removaluntreated:Area_ha:log10(tot.arth.bm) 0.568004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.734 on 482 degrees of freedom
## Multiple R-squared: 0.4486, Adjusted R-squared: 0.3983
## F-statistic: 8.913 on 44 and 482 DF, p-value: < 2.2e-16########make plots to assess the distribution of the residuals
resid_panel(bird.model.2, plots = c('resid', 'qq', 'lev', 'hist'))## Warning in plot_lev(model = model, type = type, smoother = smoother, theme = theme, : Observations with a leverage value of 1 are not included
## in the residuals versus leverage plot.
#Report F tests for terms of the model.
Anova(bird.model.2)## Note: model has aliased coefficients
## sums of squares computed by model comparison## Anova Table (Type II tests)
##
## Response: foraging.ht.m
## Sum Sq Df F value Pr(>F)
## SPECIES 158.7 5 2.2759 0.0460655
## Rat_Removal 71.1 1 5.1016 0.0243495
## Area_ha 1488.9 1 106.7906 < 2.2e-16
## log10(tot.arth.bm) 175.9 1 12.6154 0.0004202
## SPECIES:Rat_Removal 42.5 5 0.6103 0.6920877
## SPECIES:Area_ha 194.0 5 2.7834 0.0171925
## Rat_Removal:Area_ha 11.1 1 0.7934 0.3735222
## SPECIES:log10(tot.arth.bm) 36.7 5 0.5267 0.7561506
## Rat_Removal:log10(tot.arth.bm) 3.4 1 0.2473 0.6192343
## Area_ha:log10(tot.arth.bm) 1.8 1 0.1305 0.7180548
## SPECIES:Rat_Removal:Area_ha 156.8 4 2.8115 0.0250357
## SPECIES:Rat_Removal:log10(tot.arth.bm) 62.2 4 1.1150 0.3486611
## SPECIES:Area_ha:log10(tot.arth.bm) 106.4 4 1.9087 0.1077843
## Rat_Removal:Area_ha:log10(tot.arth.bm) 1133.3 1 81.2828 < 2.2e-16
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm) 182.7 4 3.2756 0.0115017
## Residuals 6720.1 482
##
## SPECIES *
## Rat_Removal *
## Area_ha ***
## log10(tot.arth.bm) ***
## SPECIES:Rat_Removal
## SPECIES:Area_ha *
## Rat_Removal:Area_ha
## SPECIES:log10(tot.arth.bm)
## Rat_Removal:log10(tot.arth.bm)
## Area_ha:log10(tot.arth.bm)
## SPECIES:Rat_Removal:Area_ha *
## SPECIES:Rat_Removal:log10(tot.arth.bm)
## SPECIES:Area_ha:log10(tot.arth.bm)
## Rat_Removal:Area_ha:log10(tot.arth.bm) ***
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm) *
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Create an effects plot displaying fitted effects
library(ggeffects)
plot(ggeffect(bird.model.2, terms = c('SPECIES', 'Rat_Removal')))## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
#How do you interpret the results at this stage?
### Taking into account Kipuku area and arthropod biomass, treatment serves as a weaker predictor for foraging height. This can be seen in the smaller F values for treatment and in the effects model. (5) One thing we have not accounted for in our models is the fact that the study design involves repeated measurements of the same kipuka (and potentially the same birds) over time. We will learn how to properly account for repeated measures, by adding a random effect for kipuka, later in the course. However, there is a different additional predictor that may be helpful in testing the treatment effect. Add DATE as a factor to the linear model from #4. How does the inclusion of the date of sampling alter the model results? Why do you think that is?
#Make a new model including date
bird$DATE <- as.factor(bird$DATE)
bird.model.3 <- lm(foraging.ht.m ~ SPECIES * Rat_Removal * Area_ha * log10(tot.arth.bm)* DATE, data = bird)
########make plots to assess the distribution of the residuals
resid_panel(bird.model.3, plots = c('resid', 'qq', 'lev', 'hist'))## Warning in plot_lev(model = model, type = type, smoother = smoother, theme = theme, : Observations with a leverage value of 1 are not included
## in the residuals versus leverage plot.
#Report F tests for terms of the model.
Anova(bird.model.3)## Note: model has aliased coefficients
## sums of squares computed by model comparison## Anova Table (Type II tests)
##
## Response: foraging.ht.m
## Sum Sq Df F value
## SPECIES 59.29 5 2.4567
## Rat_Removal 0.51 1 0.1063
## Area_ha 109.33 1 22.6487
## log10(tot.arth.bm) 77.54 1 16.0641
## DATE 2802.51 26 22.3300
## SPECIES:Rat_Removal 0
## SPECIES:Area_ha 0.10 2 0.0108
## Rat_Removal:Area_ha 1.64 1 0.3394
## SPECIES:log10(tot.arth.bm) 31.95 5 1.3236
## Rat_Removal:log10(tot.arth.bm) 0.08 1 0.0175
## Area_ha:log10(tot.arth.bm) 0.55 1 0.1136
## SPECIES:DATE 298.16 53 1.1654
## Rat_Removal:DATE 21.57 2 2.2345
## Area_ha:DATE 101.88 5 4.2210
## log10(tot.arth.bm):DATE 520.91 22 4.9051
## SPECIES:Rat_Removal:Area_ha 0
## SPECIES:Rat_Removal:log10(tot.arth.bm) 0
## SPECIES:Area_ha:log10(tot.arth.bm) 0.20 1 0.0404
## Rat_Removal:Area_ha:log10(tot.arth.bm) 9.03 1 1.8712
## SPECIES:Rat_Removal:DATE 0
## SPECIES:Area_ha:DATE 0
## Rat_Removal:Area_ha:DATE 3.27 1 0.6767
## SPECIES:log10(tot.arth.bm):DATE 31.09 13 0.4954
## Rat_Removal:log10(tot.arth.bm):DATE 53.76 1 11.1374
## Area_ha:log10(tot.arth.bm):DATE 81.47 4 4.2194
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm) 0
## SPECIES:Rat_Removal:Area_ha:DATE 0
## SPECIES:Rat_Removal:log10(tot.arth.bm):DATE 0
## SPECIES:Area_ha:log10(tot.arth.bm):DATE 0
## Rat_Removal:Area_ha:log10(tot.arth.bm):DATE 0
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm):DATE 0
## Residuals 1438.47 298
## Pr(>F)
## SPECIES 0.0334776 *
## Rat_Removal 0.7445729
## Area_ha 3.037e-06 ***
## log10(tot.arth.bm) 7.743e-05 ***
## DATE < 2.2e-16 ***
## SPECIES:Rat_Removal
## SPECIES:Area_ha 0.9892378
## Rat_Removal:Area_ha 0.5606129
## SPECIES:log10(tot.arth.bm) 0.2539311
## Rat_Removal:log10(tot.arth.bm) 0.8947990
## Area_ha:log10(tot.arth.bm) 0.7362835
## SPECIES:DATE 0.2157182
## Rat_Removal:DATE 0.1088396
## Area_ha:DATE 0.0010095 **
## log10(tot.arth.bm):DATE 4.603e-11 ***
## SPECIES:Rat_Removal:Area_ha
## SPECIES:Rat_Removal:log10(tot.arth.bm)
## SPECIES:Area_ha:log10(tot.arth.bm) 0.8407532
## Rat_Removal:Area_ha:log10(tot.arth.bm) 0.1723686
## SPECIES:Rat_Removal:DATE
## SPECIES:Area_ha:DATE
## Rat_Removal:Area_ha:DATE 0.4113965
## SPECIES:log10(tot.arth.bm):DATE 0.9265680
## Rat_Removal:log10(tot.arth.bm):DATE 0.0009533 ***
## Area_ha:log10(tot.arth.bm):DATE 0.0024472 **
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm)
## SPECIES:Rat_Removal:Area_ha:DATE
## SPECIES:Rat_Removal:log10(tot.arth.bm):DATE
## SPECIES:Area_ha:log10(tot.arth.bm):DATE
## Rat_Removal:Area_ha:log10(tot.arth.bm):DATE
## SPECIES:Rat_Removal:Area_ha:log10(tot.arth.bm):DATE
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#How does the inclusion of the date of sampling alter the model results? Why do you think that is?
### Taking into account Date as a a factor, as well as Species, treatment, area of kipuku, and arthropod biomass, makes for a pretty complicated model. The treatment effect is not significant (rat removal) when this is included in the model. Maybe this is because of the non-independence of samples, because birds could be sampled more than once at the same place at different times. Which is why we would need to do a repeated measures anova?