Data analysis for:
Foraging actively can be advantageous in heterogeneous environments
Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, United States.
Library
library(AICcmodavg)
library(emo)
library(icons)
library(magick)
library(knitr)
library(MuMIn)
library(nlme)
library(png)
library(raster)
library(scales)
library(vioplot)
library(xtable)
Loading data
data <- read.csv("data.csv")
head(data)> initial_mass final_mass strain env adult comments
> 1 0.0007245 NA s p NA died
> 2 0.0007720 NA s p NA died
> 3 0.0007095 NA s u NA died
> 4 0.0005566 NA s u NA died
> 5 0.0008130 NA s u NA died
> 6 0.0004055 0.0008495 s p NA <NA>data$growth <- data$final_mass - data$initial_mass
names(data)> [1] "initial_mass" "final_mass" "strain" "env" "adult" "comments" "growth"data <- data[ , c(3, 4, 7)]
data <- data[!is.na(data$growth), ]
data <- data[data$growth > 0, ]
table(data$strain, data$env)>
> p u
> r 29 29
> s 32 22data$strain[data$strain == "s"] <- "Sitter"
data$strain[data$strain == "r"] <- "Rover"
data$env[data$env == "p"] <- "Patchy"
data$env[data$env == "u"] <- "Uniform"
data <- data[data$growth < 0.002, ] # Removing abnormal data
## Movement data
vid.dat <- read.csv('./OneDrive-2024-08-19/Results/Results_by_ind.csv', row.names = 1, sep = ';')
str(vid.dat)> 'data.frame': 92 obs. of 25 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.976 0.949 0.919 0.993 0.945 ...
> $ Average_Speed : num 0.518 2.36 5.962 5.412 2.979 ...
> $ Average_Speed_Moving : num 0.531 2.488 6.486 5.448 3.152 ...
> $ Traveled_Dist : num 233 1062 2683 2435 1341 ...
> $ Meander : num 1423 412 1387 1112 2120 ...
> $ Traveled_Dist_Moving : num 233 1062 2683 2435 1341 ...
> $ Meander_moving : num 1423 412 1387 1112 2120 ...
> $ Exploration_absolute_value : num 65.2 1148.3 150.6 1059.7 504.2 ...
> $ Exploration_relative_value : num 0.00671 0.12435 0.01568 0.11397 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...vid.dat <- vid.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat[1:8, 1:5]> Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0 General Ind0 0 449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0 General Ind0 0 449.97mov.dat <- read.csv('movement.csv')
str(mov.dat)> 'data.frame': 92 obs. of 6 variables:
> $ mass : num 0.001991 0.001548 0.000821 0.001388 0.001131 ...
> $ strain : chr "r" "r" "s" "s" ...
> $ env : chr "p" "u" "p" "u" ...
> $ dis_trav: logi NA NA NA NA NA NA ...
> $ video : int 1 1 1 1 2 2 2 2 3 3 ...
> $ comments: logi NA NA NA NA NA NA ...head(mov.dat)> mass strain env dis_trav video comments
> 1 0.0019910 r p NA 1 NA
> 2 0.0015480 r u NA 1 NA
> 3 0.0008210 s p NA 1 NA
> 4 0.0013880 s u NA 1 NA
> 5 0.0011315 r p NA 2 NA
> 6 0.0014595 r u NA 2 NAs <- mov.dat[mov.dat$strain == "s", ]
r <- mov.dat[mov.dat$strain == "r", ]
dim(s)> [1] 46 6dim(r)> [1] 46 6ord.dat <- data.frame()
for(i in 1:46){
idx <- s[i, ]
idx2 <- r[i, ]
obj <- rbind(idx, idx2)
ord.dat <- rbind(ord.dat, obj)
}
ord.dat[1:8, ]> mass strain env dis_trav video comments
> 3 0.0008210 s p NA 1 NA
> 1 0.0019910 r p NA 1 NA
> 4 0.0013880 s u NA 1 NA
> 2 0.0015480 r u NA 1 NA
> 7 0.0014475 s p NA 2 NA
> 5 0.0011315 r p NA 2 NA
> 8 0.0012095 s u NA 2 NA
> 6 0.0014595 r u NA 2 NAvid.dat[1:8, 1:5]> Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0 General Ind0 0 449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0 General Ind0 0 449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0 General Ind0 0 449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0 General Ind0 0 449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0 General Ind0 0 449.97dim(ord.dat)> [1] 92 6vid.dat$mass <- ord.dat$mass
vid.dat$strain <- ord.dat$strain
vid.dat$env <- ord.dat$env
## Proportion of area visited
area.dat <- read.csv('Proportion-of-area-visited.csv', sep = ';')
area.dat <- area.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat$prop_area_visited <- area.dat[ , 2]
names(vid.dat)> [1] "Arena" "Sequence" "Individual"
> [4] "Beginning_seq" "End_seq" "Prop_time_lost"
> [7] "Smoothing_filter_window" "Smoothing_Polyorder" "Moving_threshold"
> [10] "Prop_time_moving" "Average_Speed" "Average_Speed_Moving"
> [13] "Traveled_Dist" "Meander" "Traveled_Dist_Moving"
> [16] "Meander_moving" "Exploration_absolute_value" "Exploration_relative_value"
> [19] "Exploration_method" "Exploration_area" "Exploration_aspect_param"
> [22] "Mean_nb_neighbours" "Prop_time_with_at_least_one_neighbour" "Mean_shortest_interind_distance"
> [25] "Mean_sum_interind_distances" "mass" "strain"
> [28] "env" "prop_area_visited"str(vid.dat)> 'data.frame': 92 obs. of 29 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.949 0.993 0.976 0.919 0.945 ...
> $ Average_Speed : num 2.36 5.412 0.518 5.962 2.979 ...
> $ Average_Speed_Moving : num 2.488 5.448 0.531 6.486 3.152 ...
> $ Traveled_Dist : num 1062 2435 233 2683 1341 ...
> $ Meander : num 412 1112 1423 1387 2120 ...
> $ Traveled_Dist_Moving : num 1062 2435 233 2683 1341 ...
> $ Meander_moving : num 412 1112 1423 1387 2120 ...
> $ Exploration_absolute_value : num 1148.3 1059.7 65.2 150.6 504.2 ...
> $ Exploration_relative_value : num 0.12435 0.11397 0.00671 0.01568 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...
> $ mass : num 0.000821 0.001991 0.001388 0.001548 0.001448 ...
> $ strain : chr "s" "r" "s" "r" ...
> $ env : chr "p" "p" "u" "u" ...
> $ prop_area_visited : num 0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...data.frame(vid.dat[1:12, 10:11])> Prop_time_moving Average_Speed
> Video SP - Apr 14 2024, 12 14 18.avi(3) 0.9486629 2.3603480
> Video RP - Apr 14 2024, 12 14 18.avi(4) 0.9934810 5.4120290
> Video SU - Apr 14 2024, 12 14 18.avi(1) 0.9757760 0.5177664
> Video RU - Apr 14 2024, 12 14 18.avi(2) 0.9191051 5.9617525
> Video SP - Apr 14 2024, 14 03 52.avi(1) 0.9452552 2.9791838
> Video RP - Apr 14 2024, 14 03 52.avi(2) 0.9984443 1.2822174
> Video SU - Apr 14 2024, 14 03 52.avi(3) 0.9691088 0.3256189
> Video RU - Apr 14 2024, 14 03 52.avi(4) 0.9758501 5.5721335
> Video SP - Apr 14 2024, 15 13 46.avi(1) 0.9748870 2.7576846
> Video RP - Apr 14 2024, 15 13 46.avi(2) 0.9659234 1.3956371
> Video SU - Apr 14 2024, 15 13 46.avi(3) 0.9698496 2.8616386
> Video RU - Apr 14 2024, 15 13 46.avi(4) 0.9392548 3.1724128
Modeling foraging behavior and
growth
##vid.dat <- vid.dat[vid.dat$Prop_time_moving > 0.7, ] ## extreme value
vid.dat$env[vid.dat$env == 'u'] <- 'Uniform'
vid.dat$env[vid.dat$env == 'p'] <- 'Patchy'
vid.dat$strain[vid.dat$strain == 'r'] <- 'Rover'
vid.dat$strain[vid.dat$strain == 's'] <- 'Sitter'
str(vid.dat)> 'data.frame': 92 obs. of 29 variables:
> $ Arena : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Sequence : chr "General" "General" "General" "General" ...
> $ Individual : chr "Ind0" "Ind0" "Ind0" "Ind0" ...
> $ Beginning_seq : num 0 0 0 0 0 0 0 0 0 0 ...
> $ End_seq : num 450 450 450 450 450 ...
> $ Prop_time_lost : num 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_filter_window : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Smoothing_Polyorder : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Moving_threshold : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Prop_time_moving : num 0.949 0.993 0.976 0.919 0.945 ...
> $ Average_Speed : num 2.36 5.412 0.518 5.962 2.979 ...
> $ Average_Speed_Moving : num 2.488 5.448 0.531 6.486 3.152 ...
> $ Traveled_Dist : num 1062 2435 233 2683 1341 ...
> $ Meander : num 412 1112 1423 1387 2120 ...
> $ Traveled_Dist_Moving : num 1062 2435 233 2683 1341 ...
> $ Meander_moving : num 412 1112 1423 1387 2120 ...
> $ Exploration_absolute_value : num 1148.3 1059.7 65.2 150.6 504.2 ...
> $ Exploration_relative_value : num 0.12435 0.11397 0.00671 0.01568 0.04999 ...
> $ Exploration_method : chr "Modern" "Modern" "Modern" "Modern" ...
> $ Exploration_area : int 10 10 10 10 10 10 10 10 10 10 ...
> $ Exploration_aspect_param : logi NA NA NA NA NA NA ...
> $ Mean_nb_neighbours : logi NA NA NA NA NA NA ...
> $ Prop_time_with_at_least_one_neighbour: logi NA NA NA NA NA NA ...
> $ Mean_shortest_interind_distance : logi NA NA NA NA NA NA ...
> $ Mean_sum_interind_distances : logi NA NA NA NA NA NA ...
> $ mass : num 0.000821 0.001991 0.001388 0.001548 0.001448 ...
> $ strain : chr "Sitter" "Rover" "Sitter" "Rover" ...
> $ env : chr "Patchy" "Patchy" "Uniform" "Uniform" ...
> $ prop_area_visited : num 0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...## Proportion of area covered
mod <- lm(prop_area_visited ~ mass*strain*env, data = vid.dat)
summary(mod)>
> Call:
> lm(formula = prop_area_visited ~ mass * strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.15443 -0.06215 -0.02083 0.06943 0.21500
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.40526 0.10699 3.788 0.000285 ***
> mass -133.89620 78.04521 -1.716 0.089917 .
> strainSitter -0.25333 0.13552 -1.869 0.065059 .
> envUniform -0.28803 0.13347 -2.158 0.033781 *
> mass:strainSitter 101.20618 101.16750 1.000 0.319999
> mass:envUniform 152.05488 95.00854 1.600 0.113257
> strainSitter:envUniform 0.04573 0.18349 0.249 0.803781
> mass:strainSitter:envUniform -6.03000 134.26488 -0.045 0.964285
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09032 on 84 degrees of freedom
> Multiple R-squared: 0.3546, Adjusted R-squared: 0.3008
> F-statistic: 6.592 on 7 and 84 DF, p-value: 3.273e-06mod1 <- lm(prop_area_visited ~ mass*strain + env, data = vid.dat)
summary(mod1)>
> Call:
> lm(formula = prop_area_visited ~ mass * strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.15481 -0.06221 -0.02851 0.07036 0.21089
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.26434 0.06370 4.150 7.73e-05 ***
> mass -34.18583 44.90920 -0.761 0.448584
> strainSitter -0.19425 0.08952 -2.170 0.032746 *
> envUniform -0.06584 0.01925 -3.420 0.000956 ***
> mass:strainSitter 71.56717 65.26893 1.096 0.275887
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09156 on 87 degrees of freedom
> Multiple R-squared: 0.3131, Adjusted R-squared: 0.2815
> F-statistic: 9.915 on 4 and 87 DF, p-value: 1.173e-06mod2 <- lm(prop_area_visited ~ strain*env, data = vid.dat)
summary(mod2)>
> Call:
> lm(formula = prop_area_visited ~ strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.13956 -0.06178 -0.03193 0.07027 0.21578
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.22457 0.01904 11.793 < 2e-16 ***
> strainSitter -0.11378 0.02693 -4.225 5.81e-05 ***
> envUniform -0.08135 0.02693 -3.021 0.0033 **
> strainSitter:envUniform 0.03087 0.03808 0.811 0.4198
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09132 on 88 degrees of freedom
> Multiple R-squared: 0.3088, Adjusted R-squared: 0.2852
> F-statistic: 13.1 on 3 and 88 DF, p-value: 3.767e-07mod3 <- lm(prop_area_visited ~ strain + env, data = vid.dat)
summary(mod3)>
> Call:
> lm(formula = prop_area_visited ~ strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.14694 -0.06400 -0.02426 0.07215 0.20806
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.21685 0.01646 13.175 < 2e-16 ***
> strainSitter -0.09835 0.01901 -5.175 1.4e-06 ***
> envUniform -0.06591 0.01901 -3.468 0.00081 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09115 on 89 degrees of freedom
> Multiple R-squared: 0.3036, Adjusted R-squared: 0.288
> F-statistic: 19.4 on 2 and 89 DF, p-value: 1.015e-07AICc(mod, mod1, mod2, mod3)> df AICc
> mod 9 -169.4912
> mod1 6 -170.9720
> mod2 5 -172.6838
> mod3 4 -174.2374vf1 <- varIdent(form = ~ 1 | strain)
vf2 <- varIdent(form = ~ 1 | env)
vf3 <- varIdent(form = ~ 1 | strain*env)
mod4 <- gls(prop_area_visited ~ strain*env, data = vid.dat, weights = vf1)
mod5 <- gls(prop_area_visited ~ strain*env, data = vid.dat, weights = vf2)
mod6 <- gls(prop_area_visited ~ strain*env, data = vid.dat, weights = vf3)
mod7 <- gls(prop_area_visited ~ strain*env, data = vid.dat)
AICc(mod4, mod5, mod6, mod7)> df AICc
> mod4 6 -146.2184
> mod5 6 -145.9677
> mod6 8 -146.7872
> mod7 5 -148.2570AIC(mod3, mod7)> Warning in AIC.default(mod3, mod7): models are not all fitted to the same number of observations> df AIC
> mod3 4 -174.6972
> mod7 5 -148.9547## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(mod3), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(mod3), col = "darkgrey", las = 1)
qqline(resid(mod3), col = "dodgerblue", lwd = 2)
## Distance traveled
mod8 <- lm(Traveled_Dist_Moving ~ mass*strain + env, data = vid.dat)
summary(mod8)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ mass * strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1009.14 -424.53 -53.04 350.55 1555.33
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 776.3 385.3 2.015 0.0470 *
> mass 362471.5 271648.8 1.334 0.1856
> strainSitter -169.6 541.5 -0.313 0.7549
> envUniform -210.1 116.5 -1.804 0.0747 .
> mass:strainSitter -287265.1 394801.6 -0.728 0.4688
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 553.8 on 87 degrees of freedom
> Multiple R-squared: 0.2534, Adjusted R-squared: 0.2191
> F-statistic: 7.381 on 4 and 87 DF, p-value: 3.629e-05mod9 <- lm(Traveled_Dist_Moving ~ strain*env, data = vid.dat)
summary(mod9)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain * env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1130.4 -411.9 -105.2 299.3 1534.0
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1201.74 115.07 10.443 <2e-16 ***
> strainSitter -436.73 162.74 -2.684 0.0087 **
> envUniform -53.19 162.74 -0.327 0.7445
> strainSitter:envUniform -278.99 230.15 -1.212 0.2287
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 551.9 on 88 degrees of freedom
> Multiple R-squared: 0.2501, Adjusted R-squared: 0.2245
> F-statistic: 9.783 on 3 and 88 DF, p-value: 1.233e-05mod10 <- lm(Traveled_Dist_Moving ~ strain + env, data = vid.dat)
summary(mod10)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1060.7 -414.0 -109.2 366.2 1603.8
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1271.49 99.92 12.725 < 2e-16 ***
> strainSitter -576.22 115.38 -4.994 2.92e-06 ***
> envUniform -192.69 115.38 -1.670 0.0984 .
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 553.3 on 89 degrees of freedom
> Multiple R-squared: 0.2376, Adjusted R-squared: 0.2204
> F-statistic: 13.87 on 2 and 89 DF, p-value: 5.727e-06mod11 <- lm(Traveled_Dist_Moving ~ strain, data = vid.dat)
summary(mod11)>
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1157.02 -495.06 -59.71 351.98 1507.44
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1175.15 82.39 14.263 < 2e-16 ***
> strainSitter -576.22 116.52 -4.945 3.51e-06 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 558.8 on 90 degrees of freedom
> Multiple R-squared: 0.2137, Adjusted R-squared: 0.2049
> F-statistic: 24.46 on 1 and 90 DF, p-value: 3.507e-06AICc(mod8, mod9, mod10, mod11)> df AICc
> mod8 6 1431.231
> mod9 5 1429.344
> mod10 4 1428.630
> mod11 3 1429.282## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(mod8), resid(mod8), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(mod8), resid(mod8), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(mod8), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(mod8), col = "darkgrey", las = 1)
qqline(resid(mod8), col = "dodgerblue", lwd = 2)
Modeling heterogeneity
## Growth rate
## png('heter.png', height = 7, width = 7, units = 'in', res = 350)
## Model selection
mod12 <- gls(growth ~ strain*env, data = data, weights = vf1)
mod13 <- gls(growth ~ strain*env, data = data, weights = vf2)
mod14 <- gls(growth ~ strain*env, data = data, weights = vf3)
mod15 <- gls(growth ~ strain*env, data = data)
AICc(mod12, mod13, mod14, mod15)> df AICc
> mod12 6 -1509.727
> mod13 6 -1520.713
> mod14 8 -1521.419
> mod15 5 -1511.646mod16 <- gls(growth ~ strain + env, data = data, weights = vf3)
mod17 <- gls(growth ~ strain, data = data, weights = vf3)
mod18 <- gls(growth ~ env, data = data, weights = vf3)
AICc(mod16, mod17, mod18)> df AICc
> mod16 7 -1540.994
> mod17 6 -1541.356
> mod18 6 -1560.452summary(mod18)> Generalized least squares fit by REML
> Model: growth ~ env
> Data: data
> AIC BIC logLik
> -1561.276 -1545.239 786.638
>
> Variance function:
> Structure: Different standard deviations per stratum
> Formula: ~1 | strain * env
> Parameter estimates:
> Sitter*Patchy Sitter*Uniform Rover*Uniform Rover*Patchy
> 1.0000000 1.3269145 1.3327763 0.6344861
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 0.0001712801 1.478526e-05 11.584516 0
> envUniform 0.0001574391 3.102320e-05 5.074884 0
>
> Correlation:
> (Intr)
> envUniform -0.477
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.80064487 -0.81033368 -0.04068966 0.74448669 2.87104053
>
> Residual standard error: 0.0001464178
> Degrees of freedom: 109 total; 107 residual## Proportion of area covered
xtable(data.frame(anova((mod3))), caption = 'Caption here...', digits = 3)> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sun Sep 15 13:56:15 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rrrrrr}
> \hline
> & Df & Sum.Sq & Mean.Sq & F.value & Pr..F. \\
> \hline
> strain & 1 & 0.222 & 0.222 & 26.778 & 0.000 \\
> env & 1 & 0.100 & 0.100 & 12.028 & 0.001 \\
> Residuals & 89 & 0.739 & 0.008 & & \\
> \hline
> \end{tabular}
> \caption{Caption here...}
> \end{table}## Growth
xtable(formatC(summary(mod18)$tTable, digits = 3, format = 'e'), caption = 'Caption here...', digits = 3)> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sun Sep 15 13:56:15 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rllll}
> \hline
> & Value & Std.Error & t-value & p-value \\
> \hline
> (Intercept) & 1.713e-04 & 1.479e-05 & 1.158e+01 & 1.334e-20 \\
> envUniform & 1.574e-04 & 3.102e-05 & 5.075e+00 & 1.644e-06 \\
> \hline
> \end{tabular}
> \caption{Caption here...}
> \end{table}
Figure 2
for(i in 1:23){
png(paste('figure', i, '.png', sep = ''), width = 7, height = 7, units = 'in', res = 350)
layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
RP <- readPNG(source = paste('imgs/RP-path/figure', i, '.png', sep = ''))
SP <- readPNG(source = paste('imgs/SP-path/figure', i, '.png', sep = ''))
RU <- readPNG(source = paste('imgs/RU-path/figure', i, '.png', sep = ''))
SU <- readPNG(source = paste('imgs/SU-path/figure', i, '.png', sep = ''))
play <- readPNG('imgs/play-button.png')
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Rovers')
par(new = TRUE)
rasterImage(RP, 0, 0, 1, 1)
mtext('A', at = 0, line = -0.5)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Sitters')
par(new = TRUE)
rasterImage(SP, 0, 0, 1, 1)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
par(new = TRUE)
rasterImage(RU, 0, 0, 1, 1)
mtext('B', at = 0, line = -0.5)
par(mar = c(2.5, 1, 1, 0))
plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
par(new = TRUE)
rasterImage(SU, 0, 0, 1, 1)
par(xpd = TRUE)
rasterImage(play, -0.08, 0, 0.02, 0.11)
dev.off()
}
Figure 1
summary(mod3)>
> Call:
> lm(formula = prop_area_visited ~ strain + env, data = vid.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.14694 -0.06400 -0.02426 0.07215 0.20806
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.21685 0.01646 13.175 < 2e-16 ***
> strainSitter -0.09835 0.01901 -5.175 1.4e-06 ***
> envUniform -0.06591 0.01901 -3.468 0.00081 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.09115 on 89 degrees of freedom
> Multiple R-squared: 0.3036, Adjusted R-squared: 0.288
> F-statistic: 19.4 on 2 and 89 DF, p-value: 1.015e-07rover.patchy <- coef(mod3)[1]
sitter.patchy <- coef(mod3)[1] + coef(mod3)[2]
rover.uniform <- coef(mod3)[1] + coef(mod3)[3]
sitter.uniform <- sitter.patchy + coef(mod3)[3]
rover.patchy> (Intercept)
> 0.2168478sitter.patchy> (Intercept)
> 0.1185rover.uniform> (Intercept)
> 0.1509348sitter.uniform> (Intercept)
> 0.05258696tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), mean)> Patchy Uniform
> Rover 0.2245652 0.14321739
> Sitter 0.1107826 0.06030435sd <- tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), sd)
sd> Patchy Uniform
> Rover 0.07829944 0.10867171
> Sitter 0.10227012 0.07042819##png("figure1-revised.png", height = 7, width = 7, units = "in", res = 360)
layout(matrix(c(0, 1, 1, 0,
0, 1, 1, 0,
0, 2, 2, 0,
0, 2, 2, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Basic boxplot
par(mgp = c(2, 1, 0))
vioplot(prop_area_visited ~ strain + env, data = vid.dat, border = NA, method = "jitter", side = "right", ylab = expression(paste("Proportion of area covered")~(mm^2)), xlab = "Environment", col = "white", las = 1, xaxt = 'n', ylim = c(0, 0.45), lineCol = 'white', rectCol = 'white')
axis(1, at = c(1.5, 3.5), labels = c('Patchy', 'Uniform'), cex = 0.8)
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
stripchart(prop_area_visited ~ strain + env, data = vid.dat, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
points(x = c(1, 2, 3, 4), y = c(rover.patchy, sitter.patchy, rover.uniform, sitter.uniform), bg = 'black', pch = 21, cex = 1.1)
segments(1, rover.patchy+sd[1, 1], 1, rover.patchy-sd[1, 1], lwd = 2.2)
segments(2, sitter.patchy+sd[2, 1], 2, sitter.patchy-sd[2, 1], lwd = 2.2)
segments(3, rover.uniform+sd[1, 2], 3, rover.uniform-sd[1, 2], lwd = 2.2)
segments(4, sitter.uniform+sd[2, 2], 4, sitter.uniform-sd[2, 2], lwd = 2.2)
legend('topleft', legend = c('Rover', 'Sitter'), cex = 1, pch = 16, col = c(alpha("purple", 0.5), alpha("orange", 0.5)), bty = 'n')
mtext("A", side = 2, at = 0.5, line = 3, las = 1, font = 1)
## Basic boxplot
par(mgp = c(2.8, 1, 0))
vioplot(growth ~ strain + env, data = data, border = NA, method = "jitter", side = "right", ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", col = "white", las = 1, xaxt = 'n', lineCol = 'white', rectCol = 'white')
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
axis(1, at = c(1.5, 3.5), labels = c('Patchy', 'Uniform'), cex = 0.8)
stripchart(growth ~ strain + env, data = data, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
means <- tapply(data$growth, list(data$env, data$strain), mean)
means> Rover Sitter
> Patchy 0.0001565889 0.0002030645
> Uniform 0.0003221207 0.0003373409points(x = c(1, 2, 3, 4), y = c(means[1, 1], means[1, 2], means[2, 1], means[2, 2]), bg = 'black', pch = 21, cex = 1.1)
segments(1, means[1, 1]+0.6344861*(mod18$sigma), 1, means[1, 1]-0.6344861*(mod18$sigma), lwd = 2.2)
segments(2, means[1, 2]+1.0000000*(mod18$sigma), 2, means[1, 2]-1.0000000*(mod18$sigma), lwd = 2.2)
segments(3, means[2, 1]+1.3327763*(mod18$sigma), 3, means[2, 1]-1.3327763*(mod18$sigma), lwd = 2.2)
segments(4, means[2, 2]+1.3269145*(mod18$sigma), 4, means[2, 2]-1.3269145*(mod18$sigma), lwd = 2.2)
legend('topleft', legend = c('Rover', 'Sitter'), cex = 1, pch = 16, col = c(alpha("purple", 0.5), alpha("orange", 0.5)), bty = 'n')
mtext("B", side = 2, at = 0.00096, line = 3, las = 1, font = 1)
##dev.off()