Skip to content

Commit

Permalink
just latitude analysis of interest
Browse files Browse the repository at this point in the history
  • Loading branch information
@rubysaltbush
rubysaltbush committed Apr 18, 2024
1 parent eccd96f commit 91449c3
Showing 1 changed file with 23 additions and 157 deletions.
180 changes: 23 additions & 157 deletions scripts/analysis/PGLS_with_latitude.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,156 +124,24 @@ pgls_allotb %>%
ylab("Floral longevity (log # days)")
ggsave("figures/symmetry_loglongevity_abs_latitude_allotb.pdf", width = 9, height = 5)

#** ttest & PGLS ----
# t-test of longevity by symmetry
ttests <- list()
ttests$allotblong <- t.test(pgls_allotb$spmean_long_days[pgls_allotb$sym_species == "1"],
pgls_allotb$spmean_long_days[pgls_allotb$sym_species == "0"])
ttests$allotblong
# Welch Two Sample t-test
#
# data: pgls_allotb$spmean_long_days[pgls_allotb$sym_species == "1"] and pgls_allotb$spmean_long_days[pgls_allotb$sym_species == "0"]
# t = 4.3647, df = 710.24, p-value = 1.462e-05
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 0.6167554 1.6252529
# sample estimates:
# mean of x mean of y
# 4.682592 3.561587

# t-test of latitude by symmetry
ttests$allotblat <- t.test(pgls_allotb$spmeanlat[pgls_allotb$sym_species == "1"],
pgls_allotb$spmeanlat[pgls_allotb$sym_species == "0"])
ttests$allotblat
# Welch Two Sample t-test
#
# data: pgls_allotb$spmeanlat[pgls_allotb$sym_species == "1"] and pgls_allotb$spmeanlat[pgls_allotb$sym_species == "0"]
# t = 1.1786, df = 921.11, p-value = 0.2389
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -0.5881533 2.3567204
# sample estimates:
# mean of x mean of y
# 26.37643 25.49215
#** PGLS ----

# run PGLS without latitude first, to check symmetry by itself
# list to store model output
pgls_models <- list()
pgls_models$allotbsymlong <- nlme::gls(log(spmean_long_days) ~
scale(as.numeric(sym_species)),
correlation = ape::corBrownian(phy = allotb,
form = ~spp),
data = pgls_allotb, method = "ML")
summary(pgls_models$allotbsymlong)
# Generalized least squares fit by maximum likelihood
# Model: log(spmean_long_days) ~ scale(as.numeric(sym_species))
# Data: pgls_allotb
# AIC BIC logLik
# 4524.62 4540.402 -2259.31
#
# Correlation Structure: corBrownian
# Formula: ~spp
# Parameter estimate(s):
# numeric(0)
#
# Coefficients:
# Value Std.Error t-value p-value
# (Intercept) 0.9563604 0.7457652 1.282388 0.1999
# scale(as.numeric(sym_species)) 0.1426238 0.0413410 3.449937 0.0006
#
# Correlation:
# (Intr)
# scale(as.numeric(sym_species)) 0.017
#
# Standardized residuals:
# Min Q1 Med Q3 Max
# -1.10113651 -0.28198308 -0.01443645 0.23146808 0.81545011
#
# Residual standard error: 3.047146
# Degrees of freedom: 1423 total; 1421 residual

# effect of (absolute) latitude on longevity?
pgls_models$allotblonglat <- nlme::gls(log(spmean_long_days) ~ scale(spmeanlatabs),
correlation = ape::corBrownian(phy = allotb,
form = ~spp),
data = pgls_allotb, method = "ML")
summary(pgls_models$allotblonglat)
# Generalized least squares fit by maximum likelihood
# Model: log(spmean_long_days) ~ scale(spmeanlatabs)
# Data: pgls_allotb
# AIC BIC logLik
# 4496.767 4512.548 -2245.383
#
# Correlation Structure: corBrownian
# Formula: ~spp
# Parameter estimate(s):
# numeric(0)
#
# Coefficients:
# Value Std.Error t-value p-value
# (Intercept) 0.8665400 0.7384333 1.173484 0.2408
# scale(spmeanlatabs) -0.2110542 0.0332768 -6.342376 0.0000
#
# Correlation:
# (Intr)
# scale(spmeanlatabs) 0.01
#
# Standardized residuals:
# Min Q1 Med Q3 Max
# -1.0204262 -0.3280829 0.0318071 0.3029884 0.8935845
#
# Residual standard error: 3.017469
# Degrees of freedom: 1423 total; 1421 residual


# yes looks like latitude has an effect on longevity
# can't do effect of latitude on symmetry bc binary variable (would need phylo logistic regression)
# but can dummy check the reverse??
pgls_models$allotbsymlat <- nlme::gls(spmeanlatabs ~ sym_species,
correlation = ape::corBrownian(phy = allotb,
form = ~spp),
data = pgls_allotb, method = "ML")
summary(pgls_models$allotbsymlat)
# Generalized least squares fit by maximum likelihood
# Model: spmeanlatabs ~ sym_species
# Data: pgls_allotb
# AIC BIC logLik
# 11242 11257.78 -5617.998
#
# Correlation Structure: corBrownian
# Formula: ~spp
# Parameter estimate(s):
# numeric(0)
#
# Coefficients:
# Value Std.Error t-value p-value
# (Intercept) 23.076481 7.901225 2.920620 0.0035
# sym_species1 -1.478082 0.940960 -1.570824 0.1164
#
# Correlation:
# (Intr)
# sym_species1 -0.021
#
# Standardized residuals:
# Min Q1 Med Q3 Max
# -0.7060763 -0.2233425 0.1016666 0.4622944 1.4485092
#
# Residual standard error: 32.28142
# Degrees of freedom: 1423 total; 1421 residual

# latitude does not vary by symmetry, which agrees with the scatter plot

# how do symmetry and latitude interact as fixed effects on longevity
# what is the effect of symmetry and latitude on longevity ?
pgls_models$allotbsymlatlong <- nlme::gls(log(spmean_long_days) ~
scale(as.numeric(sym_species))*scale(spmeanlatabs),
scale(as.numeric(sym_species))
+ scale(spmeanlatabs),
correlation = ape::corBrownian(phy = allotb,
form = ~spp),
data = pgls_allotb, method = "ML")
summary(pgls_models$allotbsymlatlong)
# Generalized least squares fit by maximum likelihood
# Model: log(spmean_long_days) ~ scale(as.numeric(sym_species)) * scale(spmeanlatabs)
# Model: log(spmean_long_days) ~ scale(as.numeric(sym_species)) + scale(spmeanlatabs)
# Data: pgls_allotb
# AIC BIC logLik
# 4458.418 4484.721 -2224.209
# 4491.906 4512.948 -2241.953
#
# Correlation Structure: corBrownian
# Formula: ~spp
Expand All @@ -282,29 +150,27 @@ summary(pgls_models$allotbsymlatlong)
#
# Coefficients:
# Value Std.Error t-value p-value
# (Intercept) 0.9536605 0.7281932 1.309626 0.1905
# scale(as.numeric(sym_species)) 0.1691634 0.0409237 4.133630 0.0000
# scale(spmeanlatabs) -0.1239595 0.0359223 -3.450767 0.0006
# scale(as.numeric(sym_species)):scale(spmeanlatabs) -0.1628957 0.0287153 -5.672776 0.0000
# (Intercept) 0.9090730 0.7370253 1.233435 0.2176
# scale(as.numeric(sym_species)) 0.1367322 0.0408635 3.346076 0.0008
# scale(spmeanlatabs) -0.2087726 0.0332832 -6.272614 0.0000
#
# Correlation:
# (Intr) sc(.(_)) scl(s)
# scale(as.numeric(sym_species)) 0.019
# scale(spmeanlatabs) 0.014 0.102
# scale(as.numeric(sym_species)):scale(spmeanlatabs) -0.011 -0.160 -0.405
#
# Standardized residuals:
# Min Q1 Med Q3 Max
# -1.115757251 -0.285272896 0.003389668 0.264009912 0.857956007
#
# Residual standard error: 2.972902
# Degrees of freedom: 1423 total; 1419 residual
# (Intr) s(.(_)
# scale(as.numeric(sym_species)) 0.017
# scale(spmeanlatabs) 0.011 0.044
#
# Standardized residuals:
# Min Q1 Med Q3 Max
# -1.007123509 -0.319472645 0.009280592 0.279001215 0.838333798
#
# Residual standard error: 3.010204
# Degrees of freedom: 1423 total; 1420 residual

# with all predictors scaled (and therefore numeric) all come out
# as significant, comparing effect sizes symmetry has strongest (+ve) effect,
# closely followed by interaction (-ve) and then latitude (-ve)
# as significant, comparing effect sizes latitude has strongest (-ve) effect,
# followed by symmetry (+ve)

rm(spp, pgls_allotb, allotb, sym_long_lat, ttests, pgls_models)
rm(spp, pgls_allotb, allotb, sym_long_lat, pgls_models)

# tried Bayesian models as per Song et al. (2022) paper, using brms package
# and running full interaction model symmetry still has significant +ve effect,
Expand Down

0 comments on commit 91449c3

Please sign in to comment.