# Clears out your memory (will erase EVERYTHING)
rm(list=ls())
#load packages ####
library(ggplot2) # Plotting package
library(dplyr) #a data management package
##
## 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
library(tidyr) # helps with data management
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ tibble 1.4.2 ✔ purrr 0.2.5
## ✔ readr 1.3.1 ✔ stringr 1.4.0
## ✔ tibble 1.4.2 ✔ forcats 0.3.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(DHARMa)
## Note that, since v0.1.6.2, DHARMa includes support for glmmTMB, but there are still a few minor limitations associatd with this package. Please see https://github.com/florianhartig/DHARMa/issues/16 for details, in particular if you use this for production.
library(glmmTMB)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:dplyr':
##
## recode
library(ggeffects)
setwd("/home/jkim19/SiteFidelity5/")
#can't find file even though location is correct. Uploaded manually instead
#import datasheet 20190613-Survey.csv, name should be "X20190613_Survey"
source("helpful_functions.R") # these are a few functions I wrote and use
data.raw <- read.csv("20190613-Survey.csv")
beetle.info <- read.csv("20190613-Beetle_info.csv")
#create base data frame, organize into less variables.
Site.fidelity <- select(data.raw, Date.Started, ScanID, X_fkPopulation_surveys, X_fkfieldID_surveys, Bracket_surveys, PartnerInteraction, Location, Date_Surveys, Time)
Site.fidelity <- rename(Site.fidelity, Pop = X_fkPopulation_surveys, FocalID = X_fkfieldID_surveys, Bracket = Bracket_surveys)
#filter in only necessary locations
levels(Site.fidelity$Location)
## [1] "" "B" "BARK" "BARK\n" "E" "I"
## [7] "I\n" "SEE NOTES" "T"
Site.fidelity <- filter(Site.fidelity, Location %sans% c("BARK", "BARK\n", "SEE NOTES", ""))
Site.fidelity <- droplevels(Site.fidelity)
#filter in only necessary bracket numbers
levels(Site.fidelity$Bracket)
## [1] "" "00"
## [3] "057" "1"
## [5] "1 " "10"
## [7] "10 " "100"
## [9] "101" "102"
## [11] "103" "104"
## [13] "105" "105 "
## [15] "106" "107"
## [17] "108" "109"
## [19] "11" "110"
## [21] "111" "112"
## [23] "113" "114"
## [25] "115" "116"
## [27] "117" "118"
## [29] "119" "12"
## [31] "120" "121"
## [33] "122" "123"
## [35] "124" "125"
## [37] "126" "127"
## [39] "128" "129"
## [41] "13" "130"
## [43] "131" "132"
## [45] "133" "134"
## [47] "135" "136"
## [49] "137" "138"
## [51] "139" "14"
## [53] "140" "141"
## [55] "142" "143"
## [57] "144" "145"
## [59] "146" "147"
## [61] "148" "149"
## [63] "15" "150"
## [65] "151" "152"
## [67] "153" "154"
## [69] "155" "156"
## [71] "157" "158"
## [73] "159" "16"
## [75] "160" "161"
## [77] "162" "163"
## [79] "164" "165"
## [81] "166" "167"
## [83] "168" "17"
## [85] "170" "171"
## [87] "172" "173"
## [89] "174" "175"
## [91] "176" "177"
## [93] "178" "179"
## [95] "18" "180"
## [97] "181" "182"
## [99] "183" "184"
## [101] "185" "185 "
## [103] "186" "187"
## [105] "188" "189"
## [107] "19" "190"
## [109] "191" "192"
## [111] "193" "194"
## [113] "195" "197"
## [115] "198" "199"
## [117] "2" "20"
## [119] "20 " "200"
## [121] "201" "202"
## [123] "203" "204"
## [125] "205" "206"
## [127] "207" "208"
## [129] "209" "21"
## [131] "210" "211"
## [133] "212" "2121"
## [135] "213" "214"
## [137] "215" "216"
## [139] "217" "218"
## [141] "219" "22"
## [143] "220" "221"
## [145] "222" "223"
## [147] "224" "225"
## [149] "226" "227"
## [151] "228" "229"
## [153] "23" "23 "
## [155] "230" "231"
## [157] "232" "233"
## [159] "234" "235"
## [161] "236" "237"
## [163] "238" "24"
## [165] "240" "241"
## [167] "242" "243"
## [169] "244" "245"
## [171] "246" "247"
## [173] "248" "249"
## [175] "25" "250"
## [177] "251" "252"
## [179] "253" "255"
## [181] "256" "258"
## [183] "259" "26"
## [185] "260" "261"
## [187] "262" "263"
## [189] "264" "265"
## [191] "266" "267"
## [193] "268" "269"
## [195] "27" "270"
## [197] "271" "272"
## [199] "273" "274"
## [201] "275" "276"
## [203] "277" "278"
## [205] "279" "28"
## [207] "280" "281"
## [209] "282" "283"
## [211] "284" "285"
## [213] "286" "287"
## [215] "288" "289"
## [217] "29" "290"
## [219] "291" "292"
## [221] "293" "294"
## [223] "295" "298"
## [225] "299" "3"
## [227] "30" "300"
## [229] "301" "302"
## [231] "303" "304"
## [233] "305" "306"
## [235] "307" "308"
## [237] "309" "31"
## [239] "310" "311"
## [241] "312" "313"
## [243] "314" "315"
## [245] "316" "317"
## [247] "318" "319"
## [249] "32" "320"
## [251] "320 " "321"
## [253] "322" "323"
## [255] "324" "325"
## [257] "326" "327"
## [259] "328" "329"
## [261] "33" "330"
## [263] "330 " "331"
## [265] "332" "333"
## [267] "334" "335"
## [269] "336" "337"
## [271] "338" "339"
## [273] "34" "340"
## [275] "341" "342"
## [277] "343" "344"
## [279] "345" "346"
## [281] "347" "348"
## [283] "349" "35"
## [285] "350" "351"
## [287] "352" "353"
## [289] "354" "355"
## [291] "356" "357"
## [293] "358" "359"
## [295] "36" "360"
## [297] "361" "362"
## [299] "363" "364"
## [301] "365" "366"
## [303] "367" "368"
## [305] "369" "37"
## [307] "370" "371"
## [309] "372" "373"
## [311] "374" "375"
## [313] "376" "377"
## [315] "378" "379"
## [317] "38" "380"
## [319] "381" "382"
## [321] "383" "384"
## [323] "385" "386"
## [325] "387" "388"
## [327] "389" "39"
## [329] "390" "391"
## [331] "392" "393"
## [333] "394" "395"
## [335] "396" "397"
## [337] "398" "399"
## [339] "4" "40"
## [341] "400" "401"
## [343] "402" "403"
## [345] "404" "405"
## [347] "406" "407"
## [349] "408" "409"
## [351] "41" "410"
## [353] "411" "412"
## [355] "413" "414"
## [357] "415" "416"
## [359] "417" "418"
## [361] "419" "42"
## [363] "420" "421"
## [365] "422" "423"
## [367] "424" "425"
## [369] "426" "427"
## [371] "428" "429"
## [373] "43" "430"
## [375] "431" "432"
## [377] "433" "434"
## [379] "435" "436"
## [381] "437" "438"
## [383] "439" "44"
## [385] "44 " "440"
## [387] "441" "442"
## [389] "444" "445"
## [391] "446" "447"
## [393] "448" "449"
## [395] "45" "450"
## [397] "451" "452"
## [399] "453" "454"
## [401] "455" "456"
## [403] "457" "458"
## [405] "459" "46"
## [407] "460" "461"
## [409] "462" "463"
## [411] "464" "465"
## [413] "466" "467"
## [415] "468" "469"
## [417] "47" "470"
## [419] "471" "472"
## [421] "473" "474"
## [423] "475" "476"
## [425] "477" "478"
## [427] "479" "48"
## [429] "480" "481"
## [431] "482" "483"
## [433] "484" "485"
## [435] "486" "487"
## [437] "488" "489"
## [439] "49" "490"
## [441] "491" "492"
## [443] "493" "494"
## [445] "495" "496"
## [447] "497" "498"
## [449] "499" "5"
## [451] "5 " "50"
## [453] "51" "52"
## [455] "53" "54"
## [457] "55" "56"
## [459] "57" "58"
## [461] "59" "6"
## [463] "60" "61"
## [465] "616" "62"
## [467] "63" "64"
## [469] "65" "66"
## [471] "67" "68"
## [473] "69" "7"
## [475] "7 " "70"
## [477] "71" "710"
## [479] "72" "73"
## [481] "74" "75"
## [483] "76" "77"
## [485] "78" "79"
## [487] "8" "80"
## [489] "81" "82"
## [491] "83" "84"
## [493] "85" "86"
## [495] "87" "88"
## [497] "89" "9"
## [499] "90" "91"
## [501] "92" "93"
## [503] "933" "94"
## [505] "95" "96"
## [507] "97" "98"
## [509] "99" "B"
## [511] "motherstick" "Motherstick"
## [513] "MOTHERSTICK" "see notes"
## [515] "See Notes" "SEE NOTES"
## [517] "T" "UK"
## [519] "unlabeled bracket in fork in log" "unlabeled E of 127"
## [521] "unlabled bracket s of 32"
Site.fidelity <- filter(Site.fidelity, Bracket %sans% c("motherstick", "Motherstick", "MOTHERSTICK", "see notes", "See Notes", "SEE NOTES", "UK", "unlabeled bracket in fork in log", "unlabeled E of 127", "unlabled bracket s of 32", "00", "", "B", "T"))
Site.fidelity <- droplevels(Site.fidelity)
#create new year column to categorize different data sets, base R
Site.fidelity$Year <- substrLeft(Site.fidelity$Date.Started, 4)
Site.fidelity <- filter(Site.fidelity, Year %in% c("2016", "2017"))
#each unique observation by eliminating partner interactiosn since those double count. Group by function keeps only the columns included
data.cooked <- group_by(Site.fidelity, Year, ScanID, Pop, FocalID, Bracket)
data.cooked <- summarize(data.cooked, Count=n())
data.cooked <- filter(data.cooked, Year %in% c("2016", "2017"))
data.cooked[data.cooked == "N/A"] <- NA
data.cooked <- na.omit(data.cooked)
#Create new column for how many times a beetle visited a unique bracket
TimesVisited <- group_by(data.cooked, Year, Pop, FocalID, Bracket)
TimesVisited <- summarize(TimesVisited, Observed_No=n())
TimesVisited$FullID <- paste(TimesVisited$Year, TimesVisited$Pop, TimesVisited$FocalID, sep = "_")
#Times visited 2
TimesVisited.2 <- group_by(data.cooked, Year, Pop, FocalID, Bracket)
TimesVisited.2 <- summarize(TimesVisited.2, Observed_No=n())
TimesVisited.2$FullID <- paste(TimesVisited.2$Year, TimesVisited.2$Pop, TimesVisited.2$FocalID, TimesVisited.2$Bracket, sep = "_")
#Create dataframe for #partner interactions for each beetle.
PartnerCount <- group_by(Site.fidelity, Year, Pop, FocalID, PartnerInteraction)
PartnerCount <- summarise(PartnerCount, No_Intercation = n())
PartnerCount<-PartnerCount[PartnerCount$PartnerInteraction!="",]
PartnerCount[PartnerCount == "N/A"] <- NA
PartnerCount <- na.omit(PartnerCount)
PartnerCount <- droplevels(PartnerCount)
#total partnerinteraction numbers for each unique beetle and bracket combination
PartnerCount2 <- group_by(PartnerCount, Year, Pop, FocalID) %>% count(wt = No_Intercation)
PartnerCount2 <- rename(PartnerCount2, No_Interaction = n)
PartnerCount2$RealID <- paste(PartnerCount2$Year, PartnerCount2$Pop, PartnerCount2$FocalID, sep = "_")
PartnerCount2 <- droplevels(PartnerCount2)
#Temp data frame groups total number of scans for each population and year combination
# VAF - Since this line below has a column named "n" it is using that as a weight when it counts and I don't think we want to do that. I think you just want a count of hte # of scans for each population and year. So you have to delete the "n" column after the first Tally before you could tally again. I added a line that is "select(-n) %>% " I also added an ungroup() for good measure.
Temp <- group_by(Site.fidelity, ScanID, Year, Pop) %>% tally() %>% select(-n) %>% ungroup() %>% group_by(Year, Pop) %>% tally()
# VAF I also delete the populations 609 and 618 for the year 2017 becuase there we only visted them once or twice. To do this I delete those lines that have less than 100 observations
Temp <- Temp %>% filter(n > 100)
#Creating Occupancy variable
SF.IO <- left_join(TimesVisited, Temp, by = c("Year", "Pop"))
SF.IO <- mutate(SF.IO, IO = (Observed_No - 1)/ (n-1))
SF.IO$FullID <- paste(SF.IO$Year, SF.IO$Pop, SF.IO$FocalID, SF.IO$Bracket, sep = "_")
#time column of each time a beetle was found on a particular bracket
tm1 <- paste(Site.fidelity$Date_Surveys, Site.fidelity$Time)
tm2 <- as.POSIXct(tm1, format = "%Y%m%d %H%M")
tm3 <- as.POSIXlt(tm1, format = "%Y%m%d %H%M")
tm.all <- data.frame(Site.fidelity, ct = tm2, lt = tm3)
#creating data frames of the first and last time a beetle was seen on a particular bracket. This will be merged to calculated the total amount of time the beetle was on the bracket.
min.br.1 <- tm.all %>% group_by(Year, Pop, Bracket, FocalID) %>% summarise(min(lt))
min.br.1$FullID <- paste(min.br.1$Year, min.br.1$Pop, min.br.1$FocalID, min.br.1$Bracket, sep = "_")
min.br.1 <- min.br.1 %>% rename("Min.time"="min(lt)")
min.br.1 <- filter(min.br.1, !is.na(Min.time))
max.br.1 <- tm.all %>% group_by(Year, Pop, Bracket, FocalID) %>% summarise(max(lt))
max.br.1$FullID <- paste(max.br.1$Year, max.br.1$Pop, max.br.1$FocalID, max.br.1$Bracket, sep = "_")
max.br.1 <- max.br.1 %>% rename("Max.time"="max(lt)")
max.br.1 <- filter(max.br.1, !is.na(Max.time))
SFi <- merge(min.br.1, max.br.1, by="FullID")
SFi<-mutate(SFi, Fi = difftime(Max.time,Min.time, units="days"))
#calculating F, the amount of time that passed since the beetle was first and last seen just in general and is not bracket specific
min.beetle.F<- tm.all %>% group_by(Year, Pop, FocalID) %>% summarise(min(lt))
min.beetle.F$FullID <- paste(min.beetle.F$Year, min.beetle.F$Pop, min.beetle.F$FocalID, sep = "_")
min.beetle.F <- min.beetle.F %>% rename("Min.time"="min(lt)")
min.beetle.F <- filter(min.beetle.F, !is.na(Min.time))
max.beetle.F<- tm.all %>% group_by(Year, Pop, FocalID) %>% summarise(max(lt))
max.beetle.F$FullID <- paste(max.beetle.F$Year, max.beetle.F$Pop, max.beetle.F$FocalID, sep = "_")
max.beetle.F <- max.beetle.F %>% rename("Max.time"="max(lt)")
max.beetle.F <- filter(max.beetle.F, !is.na(Max.time))
SF <- merge(min.beetle.F, max.beetle.F, by="FullID")
SF<-mutate(SF, Fj = difftime(Max.time,Min.time, units="days"))
#combine the dataframes SF and SFI and mutate Fi/F to calcuate permanence, the amount of time a beetle was on a particular bracket relative to the total amount of time it could have potentially been there.
SFP <- left_join(SFi, SF, by = c("Year.x", "Pop.x", "FocalID.x"))
SFP <- mutate(SFP, ITi = as.numeric(Fi)/as.numeric(Fj))
SFP <- select(SFP, Year.x, Pop.x, Bracket.x, FocalID.x, FullID.y, Fj, ITi)
SFP <- rename(SFP, Year = Year.x, Pop = Pop.x, Bracket = Bracket.x, FocalID = FocalID.x, FullID = FullID.y)
SFP$FullID <- paste(SFP$Year, SFP$Pop, SFP$FocalID, SFP$Bracket, sep = "_")
#Periodicity (It) calculation, an index to indicate the recursive likelihood of a beetle returning to a specific bracket.
SF.It <- left_join(TimesVisited, Temp, by = c("Year", "Pop"))
SF.It$FullID <- paste(SF.It$Year, SF.It$Pop, SF.It$FocalID, SF.It$Bracket, sep = "_")
SF.It <- left_join(SF.It, SFi, by = c("FullID"))
SF.It <- mutate(SF.It, It = as.numeric(Observed_No - 1)/as.numeric(Fi))
SF.It <- select(SF.It, Year, Pop, FocalID, Bracket, Observed_No, FullID, n, Fi, It)
#combine all spreadsheets including partnercount, SF.IO, SFP, SF.It variables.
##start by using FullID to merge them one by one.
Master <- left_join(SF.It, SF.IO, by = c("Year", "Pop", "FocalID", "Bracket", "Observed_No", "FullID", "n"))
Master <- merge(Master, SFP, by = "FullID")
Master <- select(Master, FullID, Year.x, Pop.x, FocalID.x, Bracket.x, Observed_No, n, Fi, Fj, It, IO, ITi)
#realid column takes out bracket from the id to combine data frames without brackets as a variable.
Master$RealID <- paste(Master$Year, Master$Pop, Master$FocalID, sep = "_")
#merging site fideilty index componenets with partner intereaction numbers
Master <- merge(Master, PartnerCount2, by = "RealID")
Master <- select(Master, RealID, FullID, Year.x, Pop.x, FocalID.x, Bracket.x, Observed_No, n, Fi, Fj, It, IO, ITi, No_Interaction)
Master <- rename(Master, Year = Year.x, Pop = Pop.x, Bracket = Bracket.x, FocalID = FocalID.x)
Master <- mutate(Master, IH1 = 3/((1/IO)+(1/ITi)+(1/It)))
str(Master)
## 'data.frame': 7531 obs. of 15 variables:
## $ RealID : chr "2016_PDR-331_HBY" "2016_PDR-331_HBY" "2016_PDR-331_HBY" "2016_PDR-331_HBY" ...
## $ FullID : chr "2016_PDR-331_HBY_199" "2016_PDR-331_HBY_6" "2016_PDR-331_HBY_11" "2016_PDR-331_HBY_13" ...
## $ Year : chr "2016" "2016" "2016" "2016" ...
## $ Pop : Factor w/ 31 levels "PDR-331","PDR-4056",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ FocalID : Factor w/ 2477 levels "","@AC","@AF",..: 264 264 264 264 265 265 265 265 265 265 ...
## $ Bracket : Factor w/ 507 levels "057","1","1 ",..: 114 460 17 39 115 496 28 90 472 448 ...
## $ Observed_No : int 1 10 4 14 1 1 5 1 3 10 ...
## $ n : int 153 153 153 153 153 153 153 153 153 153 ...
## $ Fi : 'difftime' num 0 41.3131944444444 15.3743055555556 5.99236111111111 ...
## ..- attr(*, "units")= chr "days"
## $ Fj : 'difftime' num 41.9888888888889 41.9888888888889 41.9888888888889 41.9888888888889 ...
## ..- attr(*, "units")= chr "days"
## $ It : num NaN 0.218 0.195 2.169 NaN ...
## $ IO : num 0 0.0592 0.0197 0.0855 0 ...
## $ ITi : num 0 0.984 0.366 0.143 0 ...
## $ No_Interaction: int 34 34 34 34 37 37 37 37 37 37 ...
## $ IH1 : num NaN 0.1334 0.0513 0.1566 NaN ...
#integrating beetleinfo into master spreadsheet
#beetle.info <- read.csv("/Users/johnkim/Downloads/20190613-Beetle_info.csv")
#only works ont he app, but not on the server
#import beetle info excel and simplify to mainly elytra, sex, and thoracic horn
beetle.info <- select(beetle.info, X__pkfieldID, CaptureID..InitialCap_or_Recap_or_Other, Sex, CaptureID..Population, beetle_pronotum, beetle_elytra, beetle_clypealhorn, beetle_horngap, beetle_thoracichorn, CaptureID..CaptureDate)
beetle.info$Year <- substrLeft(beetle.info$CaptureID..CaptureDate, 4)
beetle.info$TempP <- substrLeft(beetle.info$CaptureID..Population, 3)
beetle.info <- filter(beetle.info, TempP == "PDR")
beetle.info <- rename(beetle.info, Pop = CaptureID..Population, FocalID = X__pkfieldID)
beetle.info <- select(beetle.info, Year, TempP, FocalID, Sex, beetle_elytra, beetle_thoracichorn)
#merge beetleinfo into Master spreadsheet
Master <- left_join(Master, beetle.info, by = c("FocalID"))
## Warning: Column `FocalID` joining factors with different levels, coercing
## to character vector
str(Master)
## 'data.frame': 7531 obs. of 20 variables:
## $ RealID : chr "2016_PDR-331_HBY" "2016_PDR-331_HBY" "2016_PDR-331_HBY" "2016_PDR-331_HBY" ...
## $ FullID : chr "2016_PDR-331_HBY_199" "2016_PDR-331_HBY_6" "2016_PDR-331_HBY_11" "2016_PDR-331_HBY_13" ...
## $ Year.x : chr "2016" "2016" "2016" "2016" ...
## $ Pop : Factor w/ 31 levels "PDR-331","PDR-4056",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ FocalID : chr "HBY" "HBY" "HBY" "HBY" ...
## $ Bracket : Factor w/ 507 levels "057","1","1 ",..: 114 460 17 39 115 496 28 90 472 448 ...
## $ Observed_No : int 1 10 4 14 1 1 5 1 3 10 ...
## $ n : int 153 153 153 153 153 153 153 153 153 153 ...
## $ Fi : 'difftime' num 0 41.3131944444444 15.3743055555556 5.99236111111111 ...
## ..- attr(*, "units")= chr "days"
## $ Fj : 'difftime' num 41.9888888888889 41.9888888888889 41.9888888888889 41.9888888888889 ...
## ..- attr(*, "units")= chr "days"
## $ It : num NaN 0.218 0.195 2.169 NaN ...
## $ IO : num 0 0.0592 0.0197 0.0855 0 ...
## $ ITi : num 0 0.984 0.366 0.143 0 ...
## $ No_Interaction : int 34 34 34 34 37 37 37 37 37 37 ...
## $ IH1 : num NaN 0.1334 0.0513 0.1566 NaN ...
## $ Year.y : chr "2015" "2015" "2015" "2015" ...
## $ TempP : chr "PDR" "PDR" "PDR" "PDR" ...
## $ Sex : Factor w/ 3 levels "F","M","SEE NOTES": 2 2 2 2 2 2 2 2 2 2 ...
## $ beetle_elytra : num 6.66 6.66 6.66 6.66 7.18 ...
## $ beetle_thoracichorn: num 1.93 1.93 1.93 1.93 3.14 ...
Master <- select(Master, RealID, FullID, Year.x, Pop, FocalID, Bracket, Observed_No, n, Fi, Fj, It, IO, ITi, No_Interaction, IH1, Sex, beetle_elytra, beetle_thoracichorn)
Master <- rename(Master, Year = Year.x)
#calculating number times a beetle is observed for entirity of season, No_Scans gives indication for whether beetle was jsut inactive or it was just not seen much.
Site.fidelity <- group_by(Site.fidelity, Year, Pop, ScanID, FocalID) %>% tally() %>% select(-n) %>% group_by(Year, Pop, FocalID) %>% tally()
Master <- left_join(Master, Site.fidelity, by = c("Year", "Pop", "FocalID"))
## Warning: Column `FocalID` joining character vector and factor, coercing
## into character vector
Master <- rename(Master, No_Scans = n.y, n= n.x)
#must keep only 1 version of each individual beetles in form of unique Real ID by keeping each one with highest IH1 value. First convert all NA in IH1 to 0's
Master.2 <- Master
Master.2$IH1[is.na(Master.2$IH1)] <- 0
Master.2 <- Master.2 %>% group_by(RealID) %>% filter(IH1 == max(IH1))
#filter out UK and make sure no elytra is zero, check merge for sex and nonzero elytra. Filter out low No_scans to eliminate just not seen beetles. Master.3 new frame for elimintating beetles that just dont "move around" or are too idle for analysis.
Master<-Master %>% filter(FocalID %sans% c("UK", "UKF", "UKM", "UL", "ULF", "ULM" ))
Master<-Master %>% filter(!is.na(beetle_elytra))
Master.3 <- Master %>% filter(No_Scans > 2)
Master.3$IH1[is.na(Master.3$IH1)] <- 0
Master.3 <- Master.3 %>% group_by(RealID) %>% filter(IH1 == max(IH1))
Master.3 <- droplevels(Master.3)
######### statistical models ######
#graph1: IH1 vs elytra and sex
Master.3$elytra.stnd <- scale(Master.3$beetle_elytra)
Master.3$IH1.stnd <- scale(Master.3$IH1)
Master.3$No_Scans.stnd <- scale(Master.3$No_Scans)
Master.3$IH1.log_1<-log(Master.3$IH1+1)
hist(Master.3$IH1.log_1)
mod.q1 <- glmmTMB(IH1.log_1 ~ Sex * elytra.stnd + No_Scans.stnd + Year+ Pop, data = Master.3, family="gaussian")
res.mod.q1 <- simulateResiduals(mod.q1, refit = F, n= 1000)
## Model family was recoginzed or set as continous, but duplicate values were detected in the simulation - changing to integer residuals (see help for details)
newdata = Master.3
#newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod.q1, newdata = newdata))
plotResiduals(pred, res.mod.q1$scaledResiduals)
testUniformity(res.mod.q1)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.017759, p-value = 0.902
## alternative hypothesis: two-sided
testZeroInflation(res.mod.q1)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided
testDispersion(res.mod.q1)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 1.0009, p-value = 0.952
## alternative hypothesis: two.sided
Anova(mod.q1, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: IH1.log_1
## Chisq Df Pr(>Chisq)
## (Intercept) 830.3832 1 < 2.2e-16 ***
## Sex 1.4569 1 0.227420
## elytra.stnd 2.9925 1 0.083653 .
## No_Scans.stnd 1988.1440 1 < 2.2e-16 ***
## Year 1.8399 1 0.174965
## Pop 304.7032 18 < 2.2e-16 ***
## Sex:elytra.stnd 9.3819 1 0.002191 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.q1)
## Family: gaussian ( identity )
## Formula:
## IH1.log_1 ~ Sex * elytra.stnd + No_Scans.stnd + Year + Pop
## Data: Master.3
##
## AIC BIC logLik deviance df.resid
## -2590.5 -2467.1 1320.3 -2640.5 1003
##
##
## Dispersion estimate for gaussian family (sigma^2): 0.00449
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.134110 0.004654 28.82 < 2e-16 ***
## SexM 0.005197 0.004305 1.21 0.227420
## elytra.stnd -0.004976 0.002877 -1.73 0.083653 .
## No_Scans.stnd 0.097808 0.002194 44.59 < 2e-16 ***
## Year2017 0.006495 0.004789 1.36 0.174965
## PopPDR-607 -0.008694 0.007025 -1.24 0.215907
## PopPDR-609 0.002166 0.017878 0.12 0.903564
## PopPDR-611 0.062735 0.016428 3.82 0.000134 ***
## PopPDR-613 -0.040938 0.008044 -5.09 3.60e-07 ***
## PopPDR-614 0.009490 0.008721 1.09 0.276511
## PopPDR-617 0.085538 0.020674 4.14 3.51e-05 ***
## PopPDR-618 0.018083 0.015603 1.16 0.246490
## PopPDR-619 -0.014429 0.007036 -2.05 0.040295 *
## PopPDR-622 -0.040112 0.021660 -1.85 0.064048 .
## PopPDR-633 0.083425 0.013774 6.06 1.39e-09 ***
## PopPDR-670 0.003666 0.018177 0.20 0.840161
## PopPDR-674 -0.120899 0.021017 -5.75 8.80e-09 ***
## PopPDR-680 0.063152 0.009260 6.82 9.11e-12 ***
## PopPDR-681 0.051703 0.011423 4.53 6.01e-06 ***
## PopPDR-748 0.043321 0.016001 2.71 0.006780 **
## PopPDR-769 -0.098508 0.015934 -6.18 6.33e-10 ***
## PopPDR-771 0.051290 0.020912 2.45 0.014179 *
## PopPDR-795 -0.049194 0.015240 -3.23 0.001247 **
## SexM:elytra.stnd 0.013308 0.004345 3.06 0.002191 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pr.q1<-ggpredict(mod.q1, terms=c("elytra.stnd", "Sex"))
plot(pr.q1)
pr.q1$pred.unlog<-1/log(pr.q1$predicted)
pr.q1$conf.low.unlog<-1/log(pr.q1$conf.low)
pr.q1$conf.high.unlog<-1/log(pr.q1$conf.high)
pr.q1$x.unstd<-pr.q1$x+mean(Master.3$beetle_elytra)*sd(Master.3$beetle_elytra)
ggplot(pr.q1, aes(x=x.unstd, y=predicted, group=group)) +geom_line(aes(color=group))+geom_ribbon(aes(ymax=conf.high, ymin=conf.low, fill=group), alpha=0.3)+xlab( "Elytra (mm)") + ylab("Log Site Fidelity Index")+scale_color_manual(values=c('green', 'blue'), name="Sex")+scale_fill_manual(values=c('green', 'blue'), name="Sex")+cleanplot
#graph 2: Social interactions vs IH1
str(Master.3$No_Interaction)
## int [1:1028] 34 37 49 15 16 101 72 21 41 34 ...
mod.q2 <- glmmTMB(No_Interaction ~ IH1.stnd *Sex + IH1.stnd + elytra.stnd + No_Scans.stnd + Year+ Pop, data = Master.3, family="nbinom2")
res.mod.q2 <- simulateResiduals(mod.q2, refit = F, n= 1000)
newdata = Master.3
#newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod.q2, newdata = newdata))
plotResiduals(pred, res.mod.q2$scaledResiduals)
testUniformity(res.mod.q2)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.027323, p-value = 0.4266
## alternative hypothesis: two-sided
testZeroInflation(res.mod.q2)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
testDispersion(res.mod.q2)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 1.3355, p-value = 0.178
## alternative hypothesis: two.sided
Anova(mod.q2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: No_Interaction
## Chisq Df Pr(>Chisq)
## (Intercept) 4060.9397 1 < 2.2e-16 ***
## IH1.stnd 30.2134 1 3.87e-08 ***
## Sex 0.6841 1 0.408167
## elytra.stnd 0.9142 1 0.338995
## No_Scans.stnd 409.4420 1 < 2.2e-16 ***
## Year 9.8638 1 0.001686 **
## Pop 155.6679 18 < 2.2e-16 ***
## IH1.stnd:Sex 7.4480 1 0.006351 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.q2)
## Family: nbinom2 ( log )
## Formula:
## No_Interaction ~ IH1.stnd * Sex + IH1.stnd + elytra.stnd + No_Scans.stnd +
## Year + Pop
## Data: Master.3
##
## AIC BIC logLik deviance df.resid
## 6727.5 6855.8 -3337.7 6675.5 1002
##
##
## Overdispersion parameter for nbinom2 family (): 3.45
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.683312 0.042107 63.73 < 2e-16 ***
## IH1.stnd 0.218603 0.039770 5.50 3.87e-08 ***
## SexM 0.033305 0.040267 0.83 0.408167
## elytra.stnd 0.019277 0.020161 0.96 0.338995
## No_Scans.stnd 0.691481 0.034173 20.23 < 2e-16 ***
## Year2017 -0.143237 0.045607 -3.14 0.001686 **
## PopPDR-607 -0.334838 0.066116 -5.06 4.10e-07 ***
## PopPDR-609 -0.669612 0.186979 -3.58 0.000342 ***
## PopPDR-611 -0.396865 0.152592 -2.60 0.009300 **
## PopPDR-613 -0.392137 0.078386 -5.00 5.66e-07 ***
## PopPDR-614 -0.349441 0.084065 -4.16 3.23e-05 ***
## PopPDR-617 -0.439888 0.191847 -2.29 0.021853 *
## PopPDR-618 -0.637578 0.149784 -4.26 2.08e-05 ***
## PopPDR-619 -0.403817 0.066368 -6.08 1.17e-09 ***
## PopPDR-622 -0.666037 0.204323 -3.26 0.001115 **
## PopPDR-633 -0.372121 0.131124 -2.84 0.004541 **
## PopPDR-670 -0.003922 0.174261 -0.02 0.982043
## PopPDR-674 -0.685866 0.200652 -3.42 0.000630 ***
## PopPDR-680 0.159072 0.086247 1.84 0.065126 .
## PopPDR-681 0.234660 0.099785 2.35 0.018690 *
## PopPDR-748 0.008205 0.143584 0.06 0.954432
## PopPDR-769 -0.154049 0.145742 -1.06 0.290513
## PopPDR-771 -0.979463 0.227307 -4.31 1.64e-05 ***
## PopPDR-795 -0.338379 0.141172 -2.40 0.016534 *
## IH1.stnd:SexM -0.119289 0.043710 -2.73 0.006351 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pr.q2<-ggpredict(mod.q2, terms=c("IH1.stnd","Sex" ))
plot(pr.q2)
#VAF I found a call in this line that was "pr$x+mean(Master.3$IH1)*sd(Master.3$IH1" but there is no object called "pr so I changed it to pr.q2$ to the line
pr.q2$x.unstd<-pr.q2$x+mean(Master.3$IH1)*sd(Master.3$IH1)
ggplot(pr.q2, aes(x=x.unstd, y=predicted, group=group)) +geom_line(aes(color=group))+geom_ribbon(aes(ymax=conf.high, ymin=conf.low, fill=group), alpha=0.3)+xlab( "Composite Site Fidelity Index (IH1)") + ylab("Number of Social Interactions")+scale_color_manual(values=c('green', 'blue'), name="Sex")+scale_fill_manual(values=c('green', 'blue'), name="Sex")+cleanplot
#graph 3: occupancy, permanence, and periodicity
Master.4<-Master.3 %>% filter(!is.na(It))
Master.4$IO.stnd <- scale(Master.4$IO)
Master.4$It.stnd <- scale(Master.4$It)
Master.4$ITi.stnd <- scale(Master.4$ITi)
Master.4$No_Scans.stnd <- scale(Master.4$No_Scans)
Master.4$elytra.stnd <- scale(Master.4$beetle_elytra)
#graph 3: Occurrence
mod.q3 <- glmmTMB(No_Interaction ~ IO.stnd*Sex + IO.stnd + elytra.stnd + No_Scans.stnd + Year+ Pop, data = Master.4, family="nbinom2")
res.mod.q3 <- simulateResiduals(mod.q3, refit = F, n= 1000)
newdata = Master.4
#newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod.q3, newdata = newdata))
plotResiduals(pred, res.mod.q3$scaledResiduals)
testUniformity(res.mod.q3)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.018328, p-value = 0.9142
## alternative hypothesis: two-sided
testZeroInflation(res.mod.q3)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
testDispersion(res.mod.q3)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 1.183, p-value = 0.21
## alternative hypothesis: two.sided
Anova(mod.q3, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: No_Interaction
## Chisq Df Pr(>Chisq)
## (Intercept) 4436.1902 1 < 2.2e-16 ***
## IO.stnd 9.7869 1 0.001758 **
## Sex 0.1193 1 0.729822
## elytra.stnd 1.0868 1 0.297177
## No_Scans.stnd 393.6103 1 < 2.2e-16 ***
## Year 16.6189 1 4.569e-05 ***
## Pop 149.7784 18 < 2.2e-16 ***
## IO.stnd:Sex 5.6851 1 0.017109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.q3)
## Family: nbinom2 ( log )
## Formula:
## No_Interaction ~ IO.stnd * Sex + IO.stnd + elytra.stnd + No_Scans.stnd +
## Year + Pop
## Data: Master.4
##
## AIC BIC logLik deviance df.resid
## 6308.0 6433.7 -3128.0 6256.0 902
##
##
## Overdispersion parameter for nbinom2 family (): 3.69
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.84475 0.04271 66.60 < 2e-16 ***
## IO.stnd 0.12587 0.04023 3.13 0.001758 **
## SexM 0.01389 0.04023 0.35 0.729822
## elytra.stnd 0.02136 0.02049 1.04 0.297177
## No_Scans.stnd 0.67347 0.03395 19.84 < 2e-16 ***
## Year2017 -0.18588 0.04560 -4.08 4.57e-05 ***
## PopPDR-607 -0.34405 0.06567 -5.24 1.62e-07 ***
## PopPDR-609 -0.83727 0.21791 -3.84 0.000122 ***
## PopPDR-611 -0.32085 0.15717 -2.04 0.041211 *
## PopPDR-613 -0.31079 0.08362 -3.72 0.000202 ***
## PopPDR-614 -0.38440 0.08386 -4.58 4.57e-06 ***
## PopPDR-617 -0.41315 0.18549 -2.23 0.025924 *
## PopPDR-618 -0.52660 0.15610 -3.37 0.000742 ***
## PopPDR-619 -0.34094 0.06867 -4.96 6.88e-07 ***
## PopPDR-622 -0.73414 0.19963 -3.68 0.000236 ***
## PopPDR-633 -0.30216 0.12672 -2.38 0.017104 *
## PopPDR-670 -0.06697 0.17041 -0.39 0.694315
## PopPDR-674 -0.73648 0.19578 -3.76 0.000169 ***
## PopPDR-680 0.17715 0.08384 2.11 0.034595 *
## PopPDR-681 0.24551 0.09673 2.54 0.011150 *
## PopPDR-748 0.01195 0.13853 0.09 0.931263
## PopPDR-769 -0.18076 0.14198 -1.27 0.202972
## PopPDR-771 -0.97657 0.22294 -4.38 1.18e-05 ***
## PopPDR-795 -0.33508 0.13716 -2.44 0.014568 *
## IO.stnd:SexM -0.10331 0.04333 -2.38 0.017109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pr.q3<-ggpredict(mod.q3, terms=c("IO.stnd", "Sex" ))
plot(pr.q3)
ggplot(pr.q3, aes(x=x, y=predicted, group=group)) +geom_line(aes(color=group))+geom_ribbon(aes(ymax=conf.high, ymin=conf.low, fill=group), alpha=0.3)+xlab( "Occurrence (IO)") + ylab("Number of Social Interactions")+scale_color_manual(values=c('green', 'blue'), name="Sex")+scale_fill_manual(values=c('green', 'blue'), name="Sex")+cleanplot
#graph 4: Periodicity
mod.q4 <- glmmTMB(No_Interaction ~ It.stnd*Sex + It.stnd + elytra.stnd + No_Scans.stnd + Year+ Pop, data = Master.4, family="nbinom2")
res.mod.q4 <- simulateResiduals(mod.q4, refit = F, n= 1000)
newdata = Master.4
#newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod.q4, newdata = newdata))
plotResiduals(pred, res.mod.q3$scaledResiduals)
testUniformity(res.mod.q4)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.016741, p-value = 0.9572
## alternative hypothesis: two-sided
testZeroInflation(res.mod.q4)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
testDispersion(res.mod.q4)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 1.1704, p-value = 0.224
## alternative hypothesis: two.sided
Anova(mod.q4, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: No_Interaction
## Chisq Df Pr(>Chisq)
## (Intercept) 4434.8398 1 < 2.2e-16 ***
## It.stnd 0.7289 1 0.39324
## Sex 0.0370 1 0.84754
## elytra.stnd 0.6999 1 0.40280
## No_Scans.stnd 1056.8067 1 < 2.2e-16 ***
## Year 17.3015 1 3.189e-05 ***
## Pop 157.6410 18 < 2.2e-16 ***
## It.stnd:Sex 3.1001 1 0.07829 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.q4)
## Family: nbinom2 ( log )
## Formula:
## No_Interaction ~ It.stnd * Sex + It.stnd + elytra.stnd + No_Scans.stnd +
## Year + Pop
## Data: Master.4
##
## AIC BIC logLik deviance df.resid
## 6308.8 6434.5 -3128.4 6256.8 902
##
##
## Overdispersion parameter for nbinom2 family (): 3.68
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.846398 0.042742 66.59 < 2e-16 ***
## It.stnd -0.024601 0.028814 -0.85 0.393237
## SexM 0.007718 0.040146 0.19 0.847544
## elytra.stnd 0.017045 0.020373 0.84 0.402801
## No_Scans.stnd 0.729843 0.022451 32.51 < 2e-16 ***
## Year2017 -0.190047 0.045690 -4.16 3.19e-05 ***
## PopPDR-607 -0.336642 0.065674 -5.13 2.96e-07 ***
## PopPDR-609 -0.822545 0.217714 -3.78 0.000158 ***
## PopPDR-611 -0.272541 0.155070 -1.76 0.078827 .
## PopPDR-613 -0.355306 0.083468 -4.26 2.07e-05 ***
## PopPDR-614 -0.410603 0.084349 -4.87 1.13e-06 ***
## PopPDR-617 -0.332532 0.184765 -1.80 0.071899 .
## PopPDR-618 -0.476143 0.155422 -3.06 0.002187 **
## PopPDR-619 -0.338021 0.068917 -4.90 9.36e-07 ***
## PopPDR-622 -0.777135 0.200196 -3.88 0.000104 ***
## PopPDR-633 -0.178416 0.123822 -1.44 0.149611
## PopPDR-670 -0.100941 0.171068 -0.59 0.555147
## PopPDR-674 -0.797718 0.194921 -4.09 4.27e-05 ***
## PopPDR-680 0.204387 0.081997 2.49 0.012680 *
## PopPDR-681 0.247840 0.096653 2.56 0.010341 *
## PopPDR-748 0.044972 0.137722 0.33 0.744014
## PopPDR-769 -0.269613 0.140135 -1.92 0.054360 .
## PopPDR-771 -0.913870 0.223349 -4.09 4.28e-05 ***
## PopPDR-795 -0.365766 0.136971 -2.67 0.007576 **
## It.stnd:SexM -0.075431 0.042841 -1.76 0.078288 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pr.q4<-ggpredict(mod.q4, terms=c("It.stnd", "Sex"))
plot(pr.q4)
ggplot(pr.q4, aes(x=x, y=predicted, group=group)) +geom_line(aes(color=group))+geom_ribbon(aes(ymax=conf.high, ymin=conf.low, fill=group), alpha=0.3)+xlab( "Periodicity (It)") + ylab("Number of Social Interactions")+scale_color_manual(values=c('green', 'blue'), name="Sex")+scale_fill_manual(values=c('green', 'blue'), name="Sex")+cleanplot
#graph 5: Permanence
mod.q5 <- glmmTMB(No_Interaction ~ ITi.stnd*Sex + ITi.stnd + elytra.stnd + No_Scans.stnd + Year+ Pop, data = Master.4, family="nbinom2")
res.mod.q5 <- simulateResiduals(mod.q5, refit = F, n= 1000)
newdata = Master.4
#newdata$Pop = NA
newdata <- droplevels(newdata)
pred = rank(predict(mod.q5, newdata = newdata))
plotResiduals(pred, res.mod.q5$scaledResiduals)
testUniformity(res.mod.q5)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.022871, p-value = 0.7167
## alternative hypothesis: two-sided
testZeroInflation(res.mod.q5)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
testDispersion(res.mod.q5)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 1.2038, p-value = 0.208
## alternative hypothesis: two.sided
Anova(mod.q5, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: No_Interaction
## Chisq Df Pr(>Chisq)
## (Intercept) 4403.8113 1 < 2.2e-16 ***
## ITi.stnd 8.1856 1 0.004222 **
## Sex 0.0152 1 0.901750
## elytra.stnd 0.5215 1 0.470199
## No_Scans.stnd 945.4173 1 < 2.2e-16 ***
## Year 20.3407 1 6.481e-06 ***
## Pop 149.2887 18 < 2.2e-16 ***
## ITi.stnd:Sex 0.0000 1 0.997384
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod.q5)
## Family: nbinom2 ( log )
## Formula:
## No_Interaction ~ ITi.stnd * Sex + ITi.stnd + elytra.stnd + No_Scans.stnd +
## Year + Pop
## Data: Master.4
##
## AIC BIC logLik deviance df.resid
## 6305.2 6430.9 -3126.6 6253.2 902
##
##
## Overdispersion parameter for nbinom2 family (): 3.7
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.864434 0.043164 66.36 < 2e-16 ***
## ITi.stnd 0.082796 0.028939 2.86 0.004222 **
## SexM 0.004949 0.040086 0.12 0.901750
## elytra.stnd 0.014719 0.020381 0.72 0.470199
## No_Scans.stnd 0.708903 0.023055 30.75 < 2e-16 ***
## Year2017 -0.207472 0.046002 -4.51 6.48e-06 ***
## PopPDR-607 -0.335162 0.065574 -5.11 3.20e-07 ***
## PopPDR-609 -0.885009 0.218732 -4.05 5.21e-05 ***
## PopPDR-611 -0.338886 0.156808 -2.16 0.030684 *
## PopPDR-613 -0.334643 0.082930 -4.04 5.45e-05 ***
## PopPDR-614 -0.410312 0.083929 -4.89 1.01e-06 ***
## PopPDR-617 -0.415508 0.183649 -2.26 0.023666 *
## PopPDR-618 -0.521652 0.155784 -3.35 0.000812 ***
## PopPDR-619 -0.341759 0.068624 -4.98 6.35e-07 ***
## PopPDR-622 -0.782194 0.198197 -3.95 7.93e-05 ***
## PopPDR-633 -0.279611 0.124313 -2.25 0.024497 *
## PopPDR-670 -0.068791 0.170595 -0.40 0.686769
## PopPDR-674 -0.761100 0.194583 -3.91 9.17e-05 ***
## PopPDR-680 0.160539 0.083362 1.93 0.054128 .
## PopPDR-681 0.196325 0.098327 2.00 0.045863 *
## PopPDR-748 -0.024476 0.139588 -0.18 0.860809
## PopPDR-769 -0.220325 0.139828 -1.58 0.115099
## PopPDR-771 -0.965839 0.222290 -4.34 1.39e-05 ***
## PopPDR-795 -0.354190 0.136728 -2.59 0.009585 **
## ITi.stnd:SexM -0.000134 0.040866 0.00 0.997384
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pr.q5<-ggpredict(mod.q5, terms=c("ITi.stnd", "Sex"))
plot(pr.q5)
ggplot(pr.q5, aes(x=x, y=predicted, group=group)) +geom_line(aes(color=group))+geom_ribbon(aes(ymax=conf.high, ymin=conf.low, fill=group), alpha=0.3)+xlab( "Permanence (ITi)") + ylab("Number of Social Interactions")+scale_color_manual(values=c('green', 'blue'), name="Sex")+scale_fill_manual(values=c('green', 'blue'), name="Sex")+cleanplot
