I fit a rather poor quality linear model
library(tidyverse)
library(janitor)
library(gridExtra)
library(grid)
library(broom)
library(cowplot)
mapcoords <- tribble(
~habitat_points , ~lat,~long,
"HP1",20.352675, 92.194123,
"HP2",20.355748, 92.195478,
"HP3",20.362292, 92.194481,
"HP4",20.37073, 92.193485,
"HP5",20.375844, 92.19216,
"HP6",20.372682, 92.191098,
"HP7",20.365790, 92.193132,
"HP8",20.362516, 92.193397,
"HP9",20.355460, 92.192472,
"HP10",20.373804, 92.193564)
input_txt <- "habitat_points benthos.10cm3 obs.shorebirds threats
HP1 6 48 13
HP2 6 33 4
HP3 0 17 5
HP4 5 26 5
HP5 0 3 3
HP6 0 9 3
HP7 0 11 1
HP8 0 13 3
HP9 0 12 4
HP10 0 0 9
HP1 6 27 5
HP2 0 19 2
HP3 0 11 1
HP4 5 5 1
HP5 0 2 11
HP6 0 0 15
HP7 7 9 7
HP8 4 10 6
HP9 0 7 6
HP10 5 10 0
HP1 0 33 10
HP2 4 28 4
HP3 5 17 6
HP4 8 18 3
HP5 7 11 2
HP6 3 9 3
HP7 2 16 3
HP8 0 22 6
HP9 4 15 5
HP10 3 8 9
HP1 4 29 7
HP2 7 16 9
HP3 10 23 3
HP4 3 10 11
HP5 7 0 11
HP6 3 4 14
HP7 2 7 3
HP8 0 8 2
HP9 4 11 1
HP10 3 9 4
HP1 17 13 9
HP2 6 5 12
HP3 2 4 9
HP4 5 6 2
HP5 1 0 6
HP6 6 0 11
HP7 2 3 5
HP8 0 5 4
HP9 3 3 5
HP10 4 0 5
HP1 5 48 7
HP2 15 39 5
HP3 5 21 5
HP4 3 23 5
HP5 4 13 8
HP6 5 7 11
HP7 2 18 3
HP8 4 27 11
HP9 3 17 4
HP10 0 5 2"
my_temp<-tempfile()
writeLines(input_txt,my_temp)
input_df <- read_delim(my_temp,"\t") %>% clean_names() %>%
group_by(habitat_points) %>%
mutate(time=row_number(),
ratio = benthos_10cm3/obs_shorebirds) %>% ungroup() %>% mutate_all(.funs = ~ifelse(is.infinite(.),NA,.)) %>% na.omit()
input_df_2 <- input_df
input_df_2$time <- input_df$time + 1
names(input_df_2) <- paste0("lag1_",names(input_df))
input_df_2<-rename(input_df_2,
time = lag1_time,
habitat_points = lag1_habitat_points)
input_df_3 <- input_df
input_df_3$time <- input_df$time + 2
names(input_df_3) <- paste0("lag2_",names(input_df))
input_df_3<-rename(input_df_3,
time = lag2_time,
habitat_points = lag2_habitat_points)
full_df <- left_join(input_df,mapcoords ) %>% left_join(input_df_2)%>% left_join(input_df_3)
# library(leaflet)
# m <- leaflet(data=mapcoords) %>%
# addTiles() %>% # Add default OpenStreetMap map tiles
# addMarkers(lng=~long, lat=~lat, popup=~HPCODE)
ratio_them <- function(x,y){
full_df[[paste0("ratio",x,"_",y)]] <<- full_df[[x]] / full_df[[y]]
}
# ratio_them("lag1_benthos_10cm3","lag2_benthos_10cm3")
ratio_them("lag1_obs_shorebirds","lag2_obs_shorebirds")
#
# ratio_them("lag1_threats","lag2_threats")
#
# ratio_them("lag1_ratio","lag2_ratio")
lm_bird_target <- lm(obs_shorebirds ~
lat + long +
# lag1_benthos_10cm3 +
# lag1_obs_shorebirds +
lag1_threats +
lag1_ratio
# +lag2_benthos_10cm3 +
# lag2_obs_shorebirds +
# lag2_threats + lag2_ratio +
# ratiolag1_obs_shorebirds_lag2_obs_shorebirds
,data = full_df)
(smry_lm <- summary(lm_bird_target))
(glanced_lm <- broom::glance(smry_lm))
full_df$predicted_birds_lm <- predict(lm_bird_target,newdata = full_df)
to_plot_df <- full_df %>% filter(time>1) %>% select(habitat_points,
time,
predicted_birds_lm,
obs_shorebirds) %>%
pivot_longer(cols=3:4,names_to="measurement",
values_to = "value")
p1<-ggplot(data=to_plot_df,aes(x=habitat_points,y=value,color=measurement),alpha=.5) +
geom_point()
p2<-ggplot(data=to_plot_df,aes(x=time,y=value,color=measurement),alpha=.5) +
geom_point()
(gt_glm_smry<-tableGrob(glanced_lm %>% mutate_if(is.numeric,round,digits=3)))
(gt_glm<-tableGrob(tidy(smry_lm) %>% mutate_if(is.numeric,round,digits=3)))
cowplot::plot_grid(p1,
p2,
gt_glm_smry,
gt_glm,
nrow=4,
ncol=1,
labels=c("predicting birds at habitats (overplot time)",
"predicting birds at time (overplot habitats)",
"summary of linear model"))
