R Notebook
metapopultion daynmic
gettin the data
data <- read.csv('write_up_1_data.csv')
library(ggplot2)immigration and extinction
for each island I want to calculate it immigration rate and extinction rate for that I will calclate 1 divided by the avarge length of time from extinction to imigration (imigration rate) and from imigration to extinction (extinction rate)
# Load necessary library
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Convert columns x1 to x30 to binary data
data <- data %>%
mutate(across(X1:X30, ~ ifelse(. > 0, 1, 0)))
# Function to calculate the mean length of continuous sets of 1s
mean_length_continuous_ones <- function(row ,value) {
values <- as.numeric(row)
rle_values <- rle(values)
lengths_ones <- rle_values$lengths[rle_values$values == value]
if (length(lengths_ones) > 0) {
return(mean(lengths_ones))
} else {
return(NA)
}
}
# Apply the function to each row and add the result as a new column
data1 <- data %>%
rowwise() %>%
mutate(extinction_rate = 1/ mean_length_continuous_ones(c_across(X1:X30) ,1) , colnilzation_rate = 1/ mean_length_continuous_ones(c_across(X1:X30) ,0)) %>%
ungroup()
# Display the modified data fram
meta_comunity <- select(data1 ,island, size, extinction_rate, colnilzation_rate)
meta_comunityplot
ggplot(meta_comunity, aes(x = size)) +
geom_point(aes(y = extinction_rate), color = "blue") +
geom_point(aes(y = colnilzation_rate), color = "red") +
geom_smooth(aes(y = extinction_rate), method = "lm", color = "blue") +
geom_smooth(aes(y = colnilzation_rate), method = "lm", color = "red") +
labs(x = "Size", y = "Rate", title = "Effect of Size on Extinction and Colonization Rates")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
connction between island size and birth rate
model_birth <- glm(birth.rate ~ size, data = data, family = Gamma(link = "log"))
summary(model_birth)##
## Call:
## glm(formula = birth.rate ~ size, family = Gamma(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.5428379 0.0921507 -49.30 2.85e-13 ***
## size 0.0051198 0.0004564 11.22 5.49e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03696812)
##
## Null deviance: 4.88461 on 11 degrees of freedom
## Residual deviance: 0.35678 on 10 degrees of freedom
## AIC: -91.654
##
## Number of Fisher Scoring iterations: 5model_death <- glm(death.rate ~ size, data = data, family = Gamma(link = "log"))
summary(model_death)##
## Call:
## glm(formula = death.rate ~ size, family = Gamma(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.5301909 0.0382272 -92.348 5.43e-16 ***
## size -0.0001181 0.0001893 -0.624 0.547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.006361719)
##
## Null deviance: 0.065287 on 11 degrees of freedom
## Residual deviance: 0.062767 on 10 degrees of freedom
## AIC: -108.22
##
## Number of Fisher Scoring iterations: 4ploting
ggplot(data, aes(x = size)) +
geom_point(aes(y = birth.rate), color = "blue") +
geom_point(aes(y = death.rate), color = "red") +
geom_smooth(aes(y = birth.rate), method = "lm", color = "blue") +
geom_smooth(aes(y = death.rate), method = "lm", color = "red") +
labs(x = "Size", y = "Rate", title = "Effect of Size on birth and death Rates")## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
colonization / extinction probability and patch_occupancy
now I will create 2 table one for extinction and on for immigration
# Load necessary library
library(dplyr)
# Function to create new columns D1 to D29 based on the given conditions
create_new_columns <- function(df) {
for (i in 2:30) {
df <- df %>%
mutate(!!paste0("X", i-1) := case_when(
get(paste0("X", i)) == 0 & get(paste0("X", i-1)) == 1 ~ "ext",
get(paste0("X", i)) == 1 & get(paste0("X", i-1)) == 0 ~ "col",
get(paste0("X", i)) == 0 & get(paste0("X", i-1)) == 0 ~ "non",
get(paste0("X", i)) == 1 & get(paste0("X", i-1)) == 1 ~ "stay",
TRUE ~ NA_character_
))
}
return(df)
}
# Apply the function to the df data frame
df_new <- create_new_columns(data) %>%
select(island, size, birth.rate, death.rate, starts_with("X"))
# Display the new data frame
print(df_new)## island size birth.rate death.rate X1 X2 X3 X4 X5 X6 X7 X8
## 1 1 131 0.019 0.028 non non col ext non col ext non
## 2 2 101 0.016 0.033 ext non non non col ext non col
## 3 3 410 0.066 0.031 stay stay stay stay stay ext col stay
## 4 4 346 0.054 0.027 stay stay stay stay stay stay stay stay
## 5 5 4 0.009 0.029 non non non non non non non col
## 6 6 246 0.051 0.026 stay stay stay stay stay stay ext non
## 7 7 263 0.048 0.026 non col stay stay stay stay stay stay
## 8 8 97 0.016 0.026 ext non non non non non col stay
## 9 9 100 0.020 0.030 non non non non non col stay ext
## 10 10 15 0.010 0.029 non non non non non non non non
## 11 11 121 0.022 0.030 non non non non col ext non non
## 12 12 100 0.021 0.030 non non non non non col stay ext
## X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20 X21 X22 X23
## 1 col stay ext non non col ext non col ext non non non non col
## 2 ext non non col stay ext non non non col ext col ext non non
## 3 stay stay stay stay stay stay ext col stay stay stay stay stay stay stay
## 4 stay ext non col stay stay stay stay stay stay stay stay stay stay stay
## 5 ext non non non non non non non non non non col ext non non
## 6 col stay stay stay stay stay stay ext non col stay stay stay stay stay
## 7 stay stay stay ext col stay stay stay stay stay stay ext col stay stay
## 8 ext non non non non col ext non non non col stay stay ext non
## 9 non non col ext non non non col stay ext non non non non non
## 10 non non col ext non non non non non non non non col ext non
## 11 col ext non col stay stay stay ext non non non non col ext non
## 12 non non col ext non non non col stay ext non non non non non
## X24 X25 X26 X27 X28 X29 X30
## 1 stay ext non non non col 1
## 2 col ext non col ext non 0
## 3 stay stay stay stay stay stay 1
## 4 stay stay stay stay stay stay 1
## 5 col ext non non non non 0
## 6 ext col stay stay stay stay 1
## 7 stay stay stay stay ext col 1
## 8 non non non non col ext 0
## 9 non col ext col ext non 0
## 10 non non non non non non 0
## 11 non col ext non non non 0
## 12 non col ext col ext non 0write.csv(df_new, "newdata.csv", row.names=FALSE)library(tidyr)
library(dplyr)
# Assuming df is already loaded in your environment
df_long <- df_new %>%
select(-birth.rate, -death.rate ,-X30) %>%
pivot_longer(cols = starts_with("X"),
names_to = "category",
values_to = "status")
# Display the first few rows of the transformed data frame
head(df_long)df_sum <- data %>%
select(starts_with("X")) %>%
summarise(across(everything(), ~ sum(. == 1, na.rm = TRUE))) %>%
pivot_longer(cols = everything(),
names_to = "category",
values_to = "path_ocopncy") %>% filter(category != "X30") %>% mutate(path_ocopncy = path_ocopncy /12)
# Display the new data frame
full <- merge(df_sum , df_long , by ="category" )colnization_data <- full %>% filter(status == "non" | status == "col") %>% mutate(status_index = case_when(status == "col" ~1 , TRUE~0))
colnization_dataglm
I will use glm to test which factor effect the extinction and immigration probabilty
conization_model <- glm(status_index ~ path_ocopncy + size,
data = colnization_data,
family = binomial(link = "logit"))
# Display the summary of the model
summary(conization_model)##
## Call:
## glm(formula = status_index ~ path_ocopncy + size, family = binomial(link = "logit"),
## data = colnization_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.783284 0.931118 -2.989 0.0028 **
## path_ocopncy 1.117065 1.971849 0.567 0.5710
## size 0.011448 0.002783 4.114 3.89e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 217.08 on 193 degrees of freedom
## Residual deviance: 194.55 on 191 degrees of freedom
## AIC: 200.55
##
## Number of Fisher Scoring iterations: 4extinction_data <- full %>% filter(status == "stay" | status == "ext") %>% mutate(status_index = case_when(status == "ext" ~1 , TRUE~0))
extinction_model <- glm(status_index ~ path_ocopncy + size,
data = extinction_data,
family = binomial(link = "logit"))
# Display the summary of the model
summary(extinction_model)##
## Call:
## glm(formula = status_index ~ path_ocopncy + size, family = binomial(link = "logit"),
## data = extinction_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.349474 1.556899 -0.224 0.8224
## path_ocopncy 5.980190 3.031655 1.973 0.0485 *
## size -0.015527 0.002532 -6.133 8.6e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 191.10 on 153 degrees of freedom
## Residual deviance: 115.57 on 151 degrees of freedom
## AIC: 121.57
##
## Number of Fisher Scoring iterations: 5simulation
I will use agent base simulation with Gillespie event base algorithm to see the patch occupancy dynamic
create_whgith_list <- function(Simulation){
whight <- c()
for (i in 1:nrow(Simulation)) {
if (Simulation$current[i] == 1) {
whight <- c(whight, Simulation$extinction_rate[i])
} else {
whight <- c(whight, Simulation$colonization_rate[i])
}
}
return(whight)
}
step <- function(Simulation ,time){
weights = create_whgith_list(Simulation)
d_time <- rexp(1, rate = sum(weights))
random_number <- sample(1:length(weights), 1, prob=weights)[1]
Simulation$current[random_number] <- 1- Simulation$current[random_number]
return(list(Simulation = Simulation , time = time + d_time) )
}
record <- data.frame(time = 0, sum = 8)
Simulation <- mutate(meta_comunity , current = 1 )
time = 0
for (i in 1:1000) {
sim <-step(Simulation ,time)
Simulation = sim$Simulation
time = sim$time
new_row <- data.frame(time = time, sum = sum(Simulation$current))
record <- rbind(record, new_row)
if (sum(Simulation$current) <1 | time > 300){
print("end")
break
}
}## [1] "end"stable_record <- filter(record , time >100)
ggplot(record, aes(x = time, y = sum/12)) +
geom_line(color = "red") +
labs(title = "Change of island ocupncy Over Time",
x = "Time",
y = "patch ocupncy") +
theme_minimal()+ geom_hline(yintercept=0 ) +geom_hline(yintercept= mean(stable_record$sum/12),color = "blue")
the minimal model
e <- mean(meta_comunity$extinction_rate)
i <- mean(meta_comunity$colnilzation_rate)
p <- i/e
print(c("in equlibrium patch ocupancy = " , p))## [1] "in equlibrium patch ocupancy = " "0.714706602075157"library("dplyr")
data2 <- data %>% select(-island , -size , -birth.rate ,-death.rate)
sum(data2)/(30*12)## [1] 0.4416667change in patch occupancy in the data itself
# Summing columns x1 to x30 into a new data frame
# Summing each column x1 to x30 separately
sum_each_column <- colSums(data[, paste0("X", 1:30)])
sum_df <- data.frame(index = 1:length(sum_each_column), column = names(sum_each_column), sum = sum_each_column)
# Display the sums of each column
ggplot(sum_df, aes(x = index, y = sum )) +
geom_line() +
labs(title = "Change of island ocupncy Over Time",
x = "Time",
y = "patch ocupncy") +
theme_minimal()+ geom_hline(yintercept=0) +geom_hline(yintercept= mean(sum_df$sum))