0.1 get distance values (km) for the transect

Code 
#this chunk is adapted from Aguillon et al., (2022) (link to code here: https://github.com/stepfanie-aguillon/flicker-HZ-movement-Evolution2022/blob/main/2_prep-scoring-data.R)


# load packages
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code 
library(nlstools)

'nlstools' has been loaded.

IMPORTANT NOTICE: Most nonlinear regression models and data set examples
related to predictive microbiolgy have been moved to the package 'nlsMicrobio'
Code 
library(geosphere)
Warning: package 'geosphere' was built under R version 4.3.3
Code 
library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
Code 
library(ggplot2)
library(mapview)

#load dataset from publicly available repository
samps<-read.csv("https://raw.githubusercontent.com/DevonDeRaad/grosbeak.rad.hybridtransect/refs/heads/main/data/grosbeak.sample.info.csv")
#subset to only samples passing genomic filtering protocols
samps<-samps[samps$passed.genomic.filtering == "TRUE",] 

#group dataset by site locality and historic/contemporary, summarize various aspects of each locale
sub_summary <- samps %>%
  group_by(site) %>%
  summarize(samples=n(),
            mtDNA_mean = mean(na.omit(mtDNA)),
            mtDNA_se = sd(na.omit(mtDNA))/samples,
            pheno_mean = mean(na.omit(male.total/12)),
            pheno_se = sd(na.omit(male.total/12))/samples,
            ancestry_mean = mean(na.omit(mac.mel.q)),
            ancestry_se = sd(na.omit(mac.mel.q))/samples,
            lat_mean = mean(na.omit(decimallatitude)),
            long_mean = mean(na.omit(decimallongitude))) %>% 
  arrange(site)
#view the results
head(sub_summary) #looks good - (DAD, 7 Jan 2025)
# A tibble: 6 × 10
   site samples mtDNA_mean mtDNA_se pheno_mean pheno_se ancestry_mean
  <int>   <int>      <dbl>    <dbl>      <dbl>    <dbl>         <dbl>
1     0       8        NaN       NA        NaN       NA         1.00 
2     1      19          1        0          1        0         1.00 
3     2       3          1        0          1        0         1.00 
4     3      16          1        0          1        0         0.984
5     4      10          1        0          1        0         0.994
6     5       4          1        0          1        0         0.992
# ℹ 3 more variables: ancestry_se <dbl>, lat_mean <dbl>, long_mean <dbl>
Code 
#subset to only the SD transect since we don't have mtDNA and phenotype info for the anchor points
sub_summary<-sub_summary[2:12,]

# mean longitude values
mlat = mean(sub_summary$lat_mean)

# first locality
# set up in longitude, latitude
p1=c(sub_summary$long_mean[1],mlat)

# for loop to work across all localities and calculate the distance
# the loop calculates the distance between longitude values from each site using a mean latitude value
# the loop then assigns this distance to the variable "distance"
distance<-c()
for(i in 1:nrow(sub_summary)){
  distance[i] <- distm(p1,c(sub_summary$long_mean[i],mlat))/1000
}
distance
 [1]   0.0000 186.5097 223.2239 255.3094 305.0567 331.3545 355.4688 375.3938
 [9] 405.1474 436.3640 573.5127
Code 
#map the linearized transect to make sure that the distances are reasonable after being linearized
sub_summary$mlat<-rep(mlat, times=11)

mapview(sub_summary, xcol = "long_mean", ycol = "lat_mean", col.regions = c("blue"), crs = 4269, grid = FALSE) +
  mapview(sub_summary, xcol = "long_mean", ycol = "mlat", col.regions = c("red"), crs = 4269, grid = FALSE)
Code 
#add distance into the dataframe
sub_summary<-sub_summary[,c(1:10)]
sub_summary$dist<-distance

0.2 evaluate mtDNA

Code 
#This chunk is adapted from: https://github.com/stepfanie-aguillon/flicker-HZ-movement-Evolution2022/blob/main/3_nls-clines-scoring-data.R

#define the function used to fit the sigmoidal curve
# set 'maxval' = the maximum value of the measured trait
maxval=1
rhs <- function(x, c, w) {
  maxval/(1+exp((4*(x-c))/w))
}

# for plotting bootstrap
#s <- seq(0,800,length=100)

### mtDNA
# modelling the cline
mtDNA_model <- nls(mtDNA_mean ~ rhs(dist, center, width),
                   data=sub_summary,
                   start=list(center=max(sub_summary$dist)/2,width=max(sub_summary$dist)/2),
                   control = list(maxiter = 500),
                   trace=T)
1.468468    (1.97e+00): par = (286.7563 286.7563)
0.4929855   (1.19e+00): par = (419.6206 300.5443)
0.1592660   (1.10e+00): par = (383.9506 134.823)
0.05831383  (5.75e-01): par = (381.5703 44.69066)
0.04907222  (1.81e-01): par = (386.7382 54.71503)
0.04772628  (5.81e-02): par = (385.2931 50.02)
0.04755574  (2.38e-02): par = (385.8566 49.73189)
0.04751946  (1.16e-02): par = (385.8083 48.97625)
0.04750998  (6.14e-03): par = (385.8758 48.72221)
0.04750723  (3.33e-03): par = (385.8879 48.54462)
0.04750642  (1.82e-03): par = (385.9008 48.45801)
0.04750618  (9.94e-04): par = (385.9062 48.4081)
0.04750611  (5.44e-04): par = (385.9096 48.38148)
0.04750608  (2.98e-04): par = (385.9114 48.36675)
0.04750608  (1.63e-04): par = (385.9123 48.35873)
0.04750607  (8.92e-05): par = (385.9129 48.35433)
0.04750607  (4.88e-05): par = (385.9132 48.35193)
0.04750607  (2.68e-05): par = (385.9133 48.35061)
0.04750607  (1.46e-05): par = (385.9134 48.34989)
0.04750607  (8.00e-06): par = (385.9135 48.34949)
Code 
# summarizing the output
summary(mtDNA_model)

Formula: mtDNA_mean ~ rhs(dist, center, width)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
center  385.913      3.432 112.457 1.76e-15 ***
width    48.349     10.374   4.661  0.00118 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07265 on 9 degrees of freedom

Number of iterations to convergence: 19 
Achieved convergence tolerance: 8.004e-06
Code 
summarymtdna<-summary(mtDNA_model)
coef(mtDNA_model)
   center     width 
385.91346  48.34949
Code 
# calculating a confidence interval
#this approach fails to converge (not sure why), so instead we will use the bootstrapping approach to generate confidence intervals here
CI_mtDNA_model <- confint(mtDNA_model,parm=c("center","width"))
Waiting for profiling to be done...
0.04977736  (3.35e-02): par = (48.34949)
0.04969549  (1.58e-02): par = (49.42262)
0.04967736  (7.28e-03): par = (49.9355)
0.04967354  (3.31e-03): par = (50.17288)
0.04967275  (1.50e-03): par = (50.28106)
0.04967259  (6.74e-04): par = (50.33001)
0.04967256  (3.03e-04): par = (50.35209)
0.04967255  (1.37e-04): par = (50.36204)
0.04967255  (6.15e-05): par = (50.36652)
0.04967255  (2.77e-05): par = (50.36853)
0.04967255  (1.24e-05): par = (50.36944)
0.04967255  (5.56e-06): par = (50.36984)
0.05582728  (1.59e-02): par = (52.33196)
0.05580762  (6.32e-03): par = (52.91036)
0.05580452  (2.47e-03): par = (53.14176)
0.05580405  (9.59e-04): par = (53.23241)
0.05580398  (3.71e-04): par = (53.26763)
0.05580397  (1.44e-04): par = (53.28127)
0.05580396  (5.55e-05): par = (53.28655)
0.05580396  (2.15e-05): par = (53.28859)
0.05580396  (8.27e-06): par = (53.28938)
0.06576935  (1.27e-02): par = (56.25193)
0.06575578  (3.61e-03): par = (56.78345)
0.06575469  (1.01e-03): par = (56.93533)
0.06575460  (2.80e-04): par = (56.97781)
0.06575460  (7.77e-05): par = (56.98961)
0.06575460  (2.15e-05): par = (56.99289)
0.06575460  (5.99e-06): par = (56.9938)
0.07947772  (6.16e-03): par = (60.80005)
0.07947426  (9.34e-04): par = (61.09993)
0.07947418  (1.40e-04): par = (61.14556)
0.07947418  (2.09e-05): par = (61.1524)
0.07947418  (3.12e-06): par = (61.15342)
0.09698677  (2.21e-03): par = (65.45527)
0.09698628  (5.56e-05): par = (65.58032)
0.09698628  (1.36e-06): par = (65.58347)
0.1183147   (9.21e-04): par = (70.17171)
0.1183146   (7.98e-05): par = (70.23198)
0.1183146   (6.91e-06): par = (70.22675)
0.04973267  (2.33e-02): par = (48.34949)
0.04969294  (1.09e-02): par = (47.63205)
0.04968432  (5.08e-03): par = (47.30152)
0.04968244  (2.38e-03): par = (47.1477)
0.04968202  (1.12e-03): par = (47.07579)
0.04968193  (5.24e-04): par = (47.04211)
0.04968191  (2.46e-04): par = (47.02631)
0.04968191  (1.15e-04): par = (47.0189)
0.04968191  (5.41e-05): par = (47.01542)
0.04968191  (2.54e-05): par = (47.01379)
0.04968191  (1.19e-05): par = (47.01302)
0.04968191  (5.59e-06): par = (47.01266)
0.05591688  (2.85e-03): par = (45.71716)
0.05591620  (1.40e-03): par = (45.80646)
0.05591604  (6.86e-04): par = (45.85035)
0.05591600  (3.37e-04): par = (45.8719)
0.05591599  (1.65e-04): par = (45.88248)
0.05591599  (8.11e-05): par = (45.88768)
0.05591599  (3.98e-05): par = (45.89023)
0.05591599  (1.96e-05): par = (45.89148)
0.05591599  (9.58e-06): par = (45.89209)
0.06622085  (2.87e-03): par = (44.76794)
0.06622000  (1.58e-03): par = (44.66971)
0.06621975  (8.73e-04): par = (44.61565)
0.06621967  (4.81e-04): par = (44.58588)
0.06621964  (2.65e-04): par = (44.56948)
0.06621964  (1.46e-04): par = (44.56045)
0.06621963  (8.05e-05): par = (44.55547)
0.06621963  (4.43e-05): par = (44.55273)
0.06621963  (2.44e-05): par = (44.55122)
0.06621963  (1.35e-05): par = (44.55038)
0.06621963  (7.42e-06): par = (44.54992)
0.08066184  (6.51e-03): par = (43.20443)
0.08065612  (4.37e-03): par = (42.9551)
0.08065353  (2.94e-03): par = (42.78816)
0.08065236  (1.98e-03): par = (42.6762)
0.08065183  (1.33e-03): par = (42.60102)
0.08065159  (8.96e-04): par = (42.55051)
0.08065149  (6.03e-04): par = (42.51656)
0.08065144  (4.06e-04): par = (42.49373)
0.08065142  (2.73e-04): par = (42.47837)
0.08065141  (1.84e-04): par = (42.46804)
0.08065140  (1.24e-04): par = (42.46109)
0.08065140  (8.31e-05): par = (42.45642)
0.08065140  (5.59e-05): par = (42.45327)
0.08065140  (3.76e-05): par = (42.45116)
0.08065140  (2.53e-05): par = (42.44973)
0.08065140  (1.70e-05): par = (42.44877)
0.08065140  (1.15e-05): par = (42.44813)
0.08065140  (7.71e-06): par = (42.4477)
0.09929859  (1.32e-02): par = (40.34813)
0.09926563  (1.19e-02): par = (39.77391)
0.09923850  (1.09e-02): par = (39.2582)
0.09921588  (9.99e-03): par = (38.79158)
0.09919681  (9.21e-03): par = (38.36673)
0.09918058  (8.52e-03): par = (37.97782)
0.09916666  (7.92e-03): par = (37.62014)
0.09915462  (7.38e-03): par = (37.28985)
0.09914417  (6.89e-03): par = (36.98378)
0.09913502  (6.46e-03): par = (36.69927)
0.09912700  (6.06e-03): par = (36.43409)
0.09911992  (5.70e-03): par = (36.18633)
0.09911366  (5.37e-03): par = (35.95436)
0.09910810  (5.06e-03): par = (35.73677)
0.09910315  (4.78e-03): par = (35.53232)
0.09909873  (4.52e-03): par = (35.33994)
0.09909478  (4.28e-03): par = (35.15868)
0.09909124  (4.05e-03): par = (34.9877)
0.09908806  (3.84e-03): par = (34.82624)
0.09908521  (3.65e-03): par = (34.67364)
0.09908264  (3.46e-03): par = (34.52929)
0.09908032  (3.29e-03): par = (34.39264)
0.09907823  (3.12e-03): par = (34.2632)
0.09907634  (2.97e-03): par = (34.14052)
0.09907464  (2.82e-03): par = (34.02418)
0.09907310  (2.68e-03): par = (33.91381)
0.09907170  (2.55e-03): par = (33.80905)
0.09907044  (2.43e-03): par = (33.70959)
0.09906930  (2.31e-03): par = (33.61512)
0.09906826  (2.20e-03): par = (33.52537)
0.09906732  (2.10e-03): par = (33.44009)
0.09906647  (2.00e-03): par = (33.35902)
0.09906570  (1.90e-03): par = (33.28196)
0.09906500  (1.81e-03): par = (33.20868)
0.09906437  (1.72e-03): par = (33.139)
0.09906379  (1.64e-03): par = (33.07272)
0.09906327  (1.56e-03): par = (33.00967)
0.09906279  (1.49e-03): par = (32.9497)
0.09906236  (1.42e-03): par = (32.89264)
0.09906197  (1.35e-03): par = (32.83835)
0.09906162  (1.29e-03): par = (32.78669)
0.09906130  (1.23e-03): par = (32.73753)
0.09906101  (1.17e-03): par = (32.69076)
0.09906074  (1.11e-03): par = (32.64624)
0.09906050  (1.06e-03): par = (32.60388)
0.09906029  (1.01e-03): par = (32.56356)
0.09906009  (9.62e-04): par = (32.52519)
0.09905991  (9.16e-04): par = (32.48867)
0.09905975  (8.73e-04): par = (32.45392)
0.09905960  (8.31e-04): par = (32.42084)
0.09905947  (7.92e-04): par = (32.38936)
0.09905934  (7.54e-04): par = (32.35939)
0.09905923  (7.18e-04): par = (32.33088)
0.09905913  (6.84e-04): par = (32.30374)
0.09905904  (6.51e-04): par = (32.27791)
0.09905896  (6.20e-04): par = (32.25332)
0.09905889  (5.90e-04): par = (32.22993)
0.09905882  (5.62e-04): par = (32.20766)
0.09905876  (5.35e-04): par = (32.18647)
0.09905870  (5.09e-04): par = (32.1663)
0.09905865  (4.85e-04): par = (32.14711)
0.09905861  (4.62e-04): par = (32.12885)
0.09905857  (4.40e-04): par = (32.11147)
0.09905853  (4.18e-04): par = (32.09493)
0.09905850  (3.98e-04): par = (32.07919)
0.09905846  (3.79e-04): par = (32.06421)
0.09905844  (3.61e-04): par = (32.04996)
0.09905841  (3.43e-04): par = (32.0364)
0.09905839  (3.27e-04): par = (32.02349)
0.09905837  (3.11e-04): par = (32.01121)
0.09905835  (2.96e-04): par = (31.99952)
0.09905833  (2.82e-04): par = (31.9884)
0.09905832  (2.68e-04): par = (31.97782)
0.09905830  (2.55e-04): par = (31.96775)
0.09905829  (2.43e-04): par = (31.95817)
0.09905828  (2.31e-04): par = (31.94906)
0.09905827  (2.20e-04): par = (31.94038)
0.09905826  (2.09e-04): par = (31.93213)
0.09905825  (1.99e-04): par = (31.92428)
0.09905824  (1.90e-04): par = (31.91681)
0.09905824  (1.80e-04): par = (31.90971)
0.09905823  (1.72e-04): par = (31.90294)
0.09905822  (1.63e-04): par = (31.89651)
0.09905822  (1.55e-04): par = (31.89039)
0.09905821  (1.48e-04): par = (31.88457)
0.09905821  (1.41e-04): par = (31.87902)
0.09905821  (1.34e-04): par = (31.87375)
0.09905820  (1.27e-04): par = (31.86874)
0.09905820  (1.21e-04): par = (31.86397)
0.09905820  (1.15e-04): par = (31.85943)
0.09905819  (1.10e-04): par = (31.85511)
0.09905819  (1.04e-04): par = (31.851)
0.09905819  (9.93e-05): par = (31.84709)
0.09905819  (9.45e-05): par = (31.84338)
0.09905819  (9.00e-05): par = (31.83984)
0.09905819  (8.55e-05): par = (31.83647)
0.09905818  (8.14e-05): par = (31.83327)
0.09905818  (7.74e-05): par = (31.83022)
0.09905818  (7.37e-05): par = (31.82733)
0.09905818  (7.01e-05): par = (31.82457)
0.09905818  (6.67e-05): par = (31.82195)
0.09905818  (6.35e-05): par = (31.81946)
0.09905818  (6.04e-05): par = (31.81708)
0.09905818  (5.74e-05): par = (31.81482)
0.09905818  (5.47e-05): par = (31.81268)
0.09905818  (5.20e-05): par = (31.81063)
0.09905818  (4.94e-05): par = (31.80869)
0.09905817  (4.70e-05): par = (31.80684)
0.09905817  (4.47e-05): par = (31.80508)
0.09905817  (4.26e-05): par = (31.80341)
0.09905817  (4.05e-05): par = (31.80182)
0.09905817  (3.85e-05): par = (31.8003)
0.09905817  (3.66e-05): par = (31.79886)
0.09905817  (3.49e-05): par = (31.79749)
0.09905817  (3.32e-05): par = (31.79619)
0.09905817  (3.15e-05): par = (31.79495)
0.09905817  (3.01e-05): par = (31.79377)
0.09905817  (2.85e-05): par = (31.79265)
0.09905817  (2.72e-05): par = (31.79158)
0.09905817  (2.58e-05): par = (31.79057)
0.09905817  (2.46e-05): par = (31.7896)
0.09905817  (2.34e-05): par = (31.78868)
0.09905817  (2.22e-05): par = (31.78781)
0.09905817  (2.12e-05): par = (31.78698)
0.09905817  (2.01e-05): par = (31.78619)
0.09905817  (1.91e-05): par = (31.78544)
0.09905817  (1.83e-05): par = (31.78472)
0.09905817  (1.73e-05): par = (31.78404)
0.09905817  (1.65e-05): par = (31.78339)
0.09905817  (1.57e-05): par = (31.78277)
0.09905817  (1.49e-05): par = (31.78219)
0.09905817  (1.42e-05): par = (31.78163)
0.09905817  (1.35e-05): par = (31.7811)
0.09905817  (1.28e-05): par = (31.7806)
0.09905817  (1.22e-05): par = (31.78012)
0.09905817  (1.16e-05): par = (31.77966)
0.09905817  (1.11e-05): par = (31.77923)
0.09905817  (1.05e-05): par = (31.77881)
0.09905817  (9.97e-06): par = (31.77842)
0.04872684  (2.22e-02): par = (385.9135)
0.04870298  (7.31e-05): par = (386.1319)
0.04870298  (5.12e-07): par = (386.1312)
0.05486921  (3.80e-02): par = (386.4158)
0.05478115  (4.06e-03): par = (386.0233)
0.05478015  (4.56e-04): par = (385.9819)
0.05478013  (5.14e-05): par = (385.9773)
0.05478013  (5.81e-06): par = (385.9768)
0.06600758  (9.94e-02): par = (385.839)
0.06519259  (1.94e-02): par = (384.7146)
0.06516317  (3.66e-03): par = (384.5143)
0.06516213  (6.88e-04): par = (384.4772)
0.06516210  (1.30e-04): par = (384.4702)
0.06516209  (2.44e-05): par = (384.4689)
0.06516209  (4.59e-06): par = (384.4687)
0.07545952  (1.23e-01): par = (383.0361)
0.07424458  (5.84e-04): par = (381.8932)
0.07424455  (3.34e-05): par = (381.8978)
0.07424455  (1.90e-06): par = (381.8981)
0.07748713  (4.28e-02): par = (378.106)
0.07734580  (8.80e-04): par = (378.2531)
0.07734574  (4.64e-07): par = (378.2562)
8.557168    (1.57e-01): par = (360.4059)
8.150337    (1.82e-01): par = (343.654)
7.593934    (2.10e-01): par = (323.0082)
6.983748    (1.59e-01): par = (299.297)
6.719739    (2.93e-01): par = (280.4711)
4.837504    (2.69e-01): par = (213.7366)
4.123946    (2.42e-01): par = (185.5285)
3.909561    (6.09e-02): par = (164.9004)
3.896892    (2.13e-02): par = (152.5754)
3.895357    (7.77e-03): par = (141.248)
3.895149    (3.36e-03): par = (130.1262)
3.895092    (7.23e-03): par = (116.9833)
3.687915    (5.23e-01): par = (22.72239)
2.895132    (5.21e-04): par = (-76.23866)
2.895131    (1.91e-04): par = (-87.09602)
2.895131    (7.04e-05): par = (-97.94732)
2.895131    (2.59e-05): par = (-108.7964)
2.895131    (9.53e-06): par = (-119.6448)
0.04856948  (2.86e-02): par = (385.9135)
0.04853356  (2.65e-03): par = (385.6098)
0.04853325  (2.37e-04): par = (385.638)
0.04853325  (2.12e-05): par = (385.6355)
0.04853325  (1.90e-06): par = (385.6357)
0.05303670  (2.67e-03): par = (385.2439)
0.05303637  (3.33e-04): par = (385.2754)
0.05303636  (4.16e-05): par = (385.2714)
0.05303636  (5.19e-06): par = (385.2719)
0.06116547  (8.09e-03): par = (384.8834)
0.06116209  (1.26e-03): par = (384.992)
0.06116201  (1.98e-04): par = (384.975)
0.06116201  (3.10e-05): par = (384.9777)
0.06116201  (4.85e-06): par = (384.9773)
0.07294487  (8.93e-03): par = (384.6638)
0.07294013  (1.66e-03): par = (384.8032)
0.07293997  (3.10e-04): par = (384.7773)
0.07293996  (5.80e-05): par = (384.7821)
0.07293996  (1.08e-05): par = (384.7812)
0.07293996  (2.02e-06): par = (384.7814)
0.08829243  (8.21e-03): par = (384.5737)
0.08828775  (1.76e-03): par = (384.7236)
0.08828754  (3.78e-04): par = (384.6915)
0.08828753  (8.10e-05): par = (384.6984)
0.08828753  (1.74e-05): par = (384.6969)
0.08828753  (3.73e-06): par = (384.6973)
0.1071046   (7.36e-03): par = (384.6077)
0.1071002   (1.76e-03): par = (384.7652)
0.1070999   (4.23e-04): par = (384.7275)
0.1070999   (1.01e-04): par = (384.7365)
0.1070999   (2.43e-05): par = (384.7343)
0.1070999   (5.82e-06): par = (384.7349)
Code 
CI_mtDNA_model
            2.5%     97.5%
center 378.16102 393.46008
width   18.04897  85.21139
Code 
# no need to bootstrap here
#bootstrap_CI_mtDNA_model <- nlsBoot(mtDNA_model,niter=999)
#summary(bootstrap_CI_mtDNA_model)
#hist(bootstrap_CI_mtDNA_model$coefboot[,1], breaks=50)
#hist(bootstrap_CI_mtDNA_model$coefboot[,2], breaks=50)

#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=mtDNA_mean,ymin=mtDNA_mean-mtDNA_se,ymax=mtDNA_mean+mtDNA_se),color="black",shape=16) +
  geom_smooth(data=sub_summary, aes(x=dist,y=mtDNA_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="black") +
  ylim(c(-0.01,1.3)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(face="bold",size=12), axis.text=element_text(size=10,color="black"))+
  geom_rect(aes(xmin = CI_mtDNA_model[1,1], xmax = CI_mtDNA_model[1,2], ymin = 1.1, ymax = 1.2), fill = NA, color = "black") + geom_segment(aes(x = summarymtdna$parameters[1,1], y = 1.1, yend = 1.2), color = "black", lwd = 1)

0.3 evaluate genomic ancestry

Code 
#define the function used to fit the sigmoidal curve
# set 'maxval' = the maximum value of the measured trait
maxval=1
rhs <- function(x, c, w) {
  maxval/(1+exp((4*(x-c))/w))
}

### genomic ancestry
# modelling the cline
ancestry_model <- nls(ancestry_mean ~ rhs(dist, center, width),
                   data=sub_summary,
                   start=list(center=max(sub_summary$dist)/2,width=max(sub_summary$dist)/2),
                   control = list(maxiter = 500),
                   trace=T)
1.380775    (2.16e+00): par = (286.7563 286.7563)
0.4206408   (1.24e+00): par = (417.2681 309.6459)
0.1325759   (1.28e+00): par = (385.3139 147.2667)
0.04078491  (6.44e-01): par = (381.218 51.27884)
0.03131379  (1.02e-01): par = (385.7438 62.56428)
0.03104228  (1.84e-02): par = (384.9119 60.82241)
0.03103313  (3.69e-03): par = (385.0755 60.94683)
0.03103271  (9.76e-04): par = (385.0563 60.84927)
0.03103267  (3.39e-04): par = (385.0632 60.83565)
0.03103267  (1.30e-04): par = (385.0634 60.82597)
0.03103267  (5.14e-05): par = (385.0639 60.82287)
0.03103267  (2.03e-05): par = (385.0641 60.82151)
0.03103267  (8.04e-06): par = (385.0641 60.821)
Code 
# summarizing the output
summary(ancestry_model)

Formula: ancestry_mean ~ rhs(dist, center, width)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
center  385.064      2.975 129.441 4.98e-16 ***
width    60.821      9.815   6.196  0.00016 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05872 on 9 degrees of freedom

Number of iterations to convergence: 12 
Achieved convergence tolerance: 8.044e-06
Code 
summary_ancestry<-summary(ancestry_model)
coef(ancestry_model)
  center    width 
385.0641  60.8210
Code 
# calculating a confidence interval
CI_ancestry_model <- confint(ancestry_model,parm=c("center","width"))
Waiting for profiling to be done...
0.03246024  (1.94e-02): par = (60.821)
0.03244418  (6.14e-03): par = (61.4023)
0.03244258  (1.92e-03): par = (61.58734)
0.03244242  (5.98e-04): par = (61.64529)
0.03244241  (1.86e-04): par = (61.66334)
0.03244241  (5.78e-05): par = (61.66896)
0.03244241  (1.80e-05): par = (61.67071)
0.03244241  (5.60e-06): par = (61.67125)
0.03650775  (1.08e-02): par = (62.49859)
0.03650230  (3.09e-03): par = (62.85045)
0.03650185  (8.76e-04): par = (62.9515)
0.03650182  (2.48e-04): par = (62.98021)
0.03650181  (7.01e-05): par = (62.98834)
0.03650181  (1.99e-05): par = (62.99064)
0.03650181  (5.59e-06): par = (62.99129)
0.04318020  (1.06e-02): par = (64.31597)
0.04317420  (2.60e-03): par = (64.70177)
0.04317383  (6.34e-04): par = (64.79707)
0.04317381  (1.54e-04): par = (64.82031)
0.04317381  (3.74e-05): par = (64.82596)
0.04317381  (9.06e-06): par = (64.82733)
0.05242306  (9.43e-03): par = (66.67866)
0.05241750  (1.83e-03): par = (67.06975)
0.05241729  (3.49e-04): par = (67.14569)
0.05241729  (6.65e-05): par = (67.16021)
0.05241729  (1.27e-05): par = (67.16298)
0.05241729  (2.40e-06): par = (67.16351)
0.06420314  (7.62e-03): par = (69.53126)
0.06419892  (1.00e-03): par = (69.89254)
0.06419885  (1.30e-04): par = (69.9402)
0.06419885  (1.68e-05): par = (69.94638)
0.06419885  (2.17e-06): par = (69.94717)
0.07850027  (5.79e-03): par = (72.78158)
0.07849747  (3.75e-04): par = (73.09505)
0.07849746  (2.38e-05): par = (73.11542)
0.07849746  (1.50e-06): par = (73.11671)
0.03244495  (8.84e-03): par = (60.821)
0.03244153  (3.06e-03): par = (60.55973)
0.03244112  (1.06e-03): par = (60.46962)
0.03244107  (3.69e-04): par = (60.43836)
0.03244106  (1.29e-04): par = (60.42749)
0.03244106  (4.47e-05): par = (60.42371)
0.03244106  (1.55e-05): par = (60.42239)
0.03244106  (5.42e-06): par = (60.42194)
0.03650518  (7.91e-03): par = (60.03344)
0.03650205  (2.92e-03): par = (60.27927)
0.03650163  (1.08e-03): par = (60.37029)
0.03650157  (3.95e-04): par = (60.40382)
0.03650156  (1.45e-04): par = (60.41615)
0.03650156  (5.33e-05): par = (60.42067)
0.03650156  (1.96e-05): par = (60.42234)
0.03650156  (7.18e-06): par = (60.42295)
0.04320461  (7.16e-03): par = (60.42397)
0.04320153  (2.81e-03): par = (60.66845)
0.04320105  (1.10e-03): par = (60.7648)
0.04320098  (4.31e-04): par = (60.80259)
0.04320097  (1.69e-04): par = (60.81738)
0.04320097  (6.60e-05): par = (60.82317)
0.04320096  (2.58e-05): par = (60.82543)
0.04320096  (1.01e-05): par = (60.82631)
0.04320096  (3.92e-06): par = (60.82666)
0.05252883  (7.71e-03): par = (61.23233)
0.05252441  (3.23e-03): par = (61.5281)
0.05252363  (1.35e-03): par = (61.6526)
0.05252350  (5.64e-04): par = (61.7047)
0.05252347  (2.35e-04): par = (61.72645)
0.05252347  (9.79e-05): par = (61.73551)
0.05252347  (4.08e-05): par = (61.73929)
0.05252347  (1.70e-05): par = (61.74086)
0.05252347  (7.08e-06): par = (61.74152)
0.06444487  (9.12e-03): par = (62.66262)
0.06443716  (4.02e-03): par = (63.06182)
0.06443566  (1.76e-03): par = (63.23861)
0.06443537  (7.71e-04): par = (63.31623)
0.06443532  (3.36e-04): par = (63.35018)
0.06443531  (1.47e-04): par = (63.365)
0.06443531  (6.40e-05): par = (63.37147)
0.06443531  (2.79e-05): par = (63.3743)
0.06443531  (1.22e-05): par = (63.37553)
0.06443531  (5.30e-06): par = (63.37606)
0.03200357  (1.38e-02): par = (385.0641)
0.03199742  (7.76e-05): par = (385.1827)
0.03199742  (3.53e-07): par = (385.1833)
0.03595026  (5.08e-03): par = (385.3235)
0.03594927  (3.29e-04): par = (385.2795)
0.03594927  (2.14e-05): par = (385.2766)
0.03594927  (1.40e-06): par = (385.2764)
0.04307271  (1.79e-02): par = (385.3635)
0.04305686  (2.56e-03): par = (385.2001)
0.04305654  (3.70e-04): par = (385.1768)
0.04305653  (5.34e-05): par = (385.1734)
0.04305653  (7.71e-06): par = (385.1729)
0.05321938  (3.84e-02): par = (385.0773)
0.05312241  (8.71e-03): par = (384.7035)
0.05311746  (1.97e-03): par = (384.6199)
0.05311721  (4.43e-04): par = (384.6011)
0.05311719  (1.00e-04): par = (384.5969)
0.05311719  (2.26e-05): par = (384.5959)
0.05311719  (5.08e-06): par = (384.5957)
0.06534579  (6.17e-02): par = (384.0543)
0.06503292  (1.44e-02): par = (383.4412)
0.06501641  (3.24e-03): par = (383.3039)
0.06501558  (7.28e-04): par = (383.2732)
0.06501554  (1.63e-04): par = (383.2663)
0.06501553  (3.66e-05): par = (383.2648)
0.06501553  (8.22e-06): par = (383.2644)
0.03190540  (1.31e-02): par = (385.0641)
0.03190028  (8.76e-04): par = (384.9422)
0.03190025  (5.78e-05): par = (384.9503)
0.03190025  (3.81e-06): par = (384.9498)
0.03514126  (3.01e-03): par = (384.808)
0.03514097  (2.98e-04): par = (384.8389)
0.03514097  (2.95e-05): par = (384.8358)
0.03514097  (2.92e-06): par = (384.8361)
0.04076398  (4.23e-03): par = (384.7162)
0.04076334  (5.43e-04): par = (384.7651)
0.04076333  (6.98e-05): par = (384.7589)
0.04076333  (8.96e-06): par = (384.7597)
0.04876756  (4.54e-03): par = (384.6788)
0.04876671  (7.04e-04): par = (384.7391)
0.04876669  (1.09e-04): par = (384.7297)
0.04876669  (1.70e-05): par = (384.7312)
0.04876669  (2.65e-06): par = (384.731)
0.05912512  (4.58e-03): par = (384.7006)
0.05912410  (8.23e-04): par = (384.7709)
0.05912407  (1.48e-04): par = (384.7582)
0.05912407  (2.67e-05): par = (384.7605)
0.05912407  (4.80e-06): par = (384.7601)
0.07180075  (4.49e-03): par = (384.791)
0.07179959  (9.06e-04): par = (384.8708)
0.07179955  (1.83e-04): par = (384.8547)
0.07179955  (3.70e-05): par = (384.858)
0.07179955  (7.49e-06): par = (384.8573)
Code 
CI_ancestry_model
            2.5%    97.5%
center 378.46299 391.6501
width   37.10079  90.8571
Code 
# bootstrapping the output
#bootstrap_CI_ancestry_model <- nlsBoot(ancestry_model,niter=999)
#summary(bootstrap_CI_ancestry_model)
#hist(bootstrap_CI_ancestry_model$coefboot[,1], breaks=50)
#hist(bootstrap_CI_ancestry_model$coefboot[,2], breaks=50)

#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=ancestry_mean,ymin=ancestry_mean-ancestry_se,ymax=ancestry_mean+ancestry_se),color="gray",shape=17) +
  geom_smooth(data=sub_summary, aes(x=dist,y=ancestry_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="gray") +
  ylim(c(-0.01,1.3)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(face="bold",size=12), axis.text=element_text(size=10,color="black"))+
  geom_rect(aes(xmin = CI_ancestry_model[1,1], xmax = CI_ancestry_model[1,2], ymin = 1.1, ymax = 1.2), fill = NA, color = "gray") + geom_segment(aes(x = summary_ancestry$parameters[1,1], y = 1.1, yend = 1.2), color = "gray", lwd = 1)

0.4 evaluate phenotype

Code 
#define the function used to fit the sigmoidal curve
# set 'maxval' = the maximum value of the measured trait
maxval=1
rhs <- function(x, c, w) {
  maxval/(1+exp((4*(x-c))/w))
}

### genomic ancestry
# modelling the cline
pheno_model <- nls(pheno_mean ~ rhs(dist, center, width),
                   data=sub_summary,
                   start=list(center=max(sub_summary$dist)/2,width=max(sub_summary$dist)/2),
                   control = list(maxiter = 500),
                   trace=T)
1.548181    (2.60e+00): par = (286.7563 286.7563)
0.4004653   (1.01e+00): par = (425.2626 360.0715)
0.2231170   (1.33e+00): par = (375.0472 10.83626)
0.1511405   (1.56e+00): par = (373.636 84.61465)
0.07411963  (7.84e-01): par = (382.8516 87.69749)
0.04441942  (2.35e-01): par = (391.6904 76.0001)
0.04112654  (9.75e-02): par = (390.4989 67.05909)
0.04054233  (4.79e-02): par = (390.5653 63.61649)
0.04040112  (2.45e-02): par = (390.5401 61.96988)
0.04036410  (1.28e-02): par = (390.5327 61.14442)
0.04035406  (6.70e-03): par = (390.5288 60.71879)
0.04035130  (3.53e-03): par = (390.5269 60.49654)
0.04035053  (1.86e-03): par = (390.5259 60.37973)
0.04035031  (9.85e-04): par = (390.5253 60.31813)
0.04035025  (5.21e-04): par = (390.5251 60.28559)
0.04035024  (2.76e-04): par = (390.5249 60.26839)
0.04035023  (1.46e-04): par = (390.5248 60.25929)
0.04035023  (7.72e-05): par = (390.5248 60.25447)
0.04035023  (4.08e-05): par = (390.5248 60.25192)
0.04035023  (2.16e-05): par = (390.5248 60.25057)
0.04035023  (1.15e-05): par = (390.5248 60.24986)
0.04035023  (6.08e-06): par = (390.5248 60.24948)
Code 
# summarizing the output
summary(pheno_model)

Formula: pheno_mean ~ rhs(dist, center, width)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
center  390.525      3.437 113.615 1.61e-15 ***
width    60.249     11.031   5.462    4e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06696 on 9 degrees of freedom

Number of iterations to convergence: 21 
Achieved convergence tolerance: 6.079e-06
Code 
summary_pheno<-summary(pheno_model)
coef(pheno_model)
   center     width 
390.52475  60.24948
Code 
# calculating a confidence interval
CI_pheno_model <- confint(pheno_model,parm=c("center","width"))
Waiting for profiling to be done...
0.04221480  (6.69e-03): par = (60.24948)
0.04221192  (3.48e-03): par = (60.02536)
0.04221114  (1.81e-03): par = (59.90932)
0.04221093  (9.42e-04): par = (59.84904)
0.04221088  (4.91e-04): par = (59.81767)
0.04221086  (2.56e-04): par = (59.80133)
0.04221086  (1.33e-04): par = (59.79282)
0.04221086  (6.95e-05): par = (59.78838)
0.04221086  (3.62e-05): par = (59.78607)
0.04221086  (1.89e-05): par = (59.78486)
0.04221086  (9.85e-06): par = (59.78424)
0.04760959  (1.16e-02): par = (59.3349)
0.04759978  (6.07e-03): par = (59.74502)
0.04759712  (3.15e-03): par = (59.95987)
0.04759640  (1.63e-03): par = (60.07169)
0.04759621  (8.45e-04): par = (60.12969)
0.04759616  (4.37e-04): par = (60.15971)
0.04759615  (2.26e-04): par = (60.17524)
0.04759614  (1.17e-04): par = (60.18327)
0.04759614  (6.03e-05): par = (60.18742)
0.04759614  (3.11e-05): par = (60.18957)
0.04759614  (1.61e-05): par = (60.19067)
0.04759614  (8.31e-06): par = (60.19125)
0.05646077  (1.15e-02): par = (60.59508)
0.05644959  (5.88e-03): par = (61.04399)
0.05644665  (3.00e-03): par = (61.27538)
0.05644588  (1.52e-03): par = (61.39367)
0.05644569  (7.73e-04): par = (61.45387)
0.05644564  (3.92e-04): par = (61.48444)
0.05644562  (1.99e-04): par = (61.49995)
0.05644562  (1.01e-04): par = (61.50781)
0.05644562  (5.10e-05): par = (61.51179)
0.05644562  (2.59e-05): par = (61.51381)
0.05644562  (1.31e-05): par = (61.51483)
0.05644562  (6.63e-06): par = (61.51535)
0.06869372  (1.19e-02): par = (62.83679)
0.06867927  (5.80e-03): par = (63.36939)
0.06867585  (2.80e-03): par = (63.62998)
0.06867505  (1.34e-03): par = (63.75604)
0.06867487  (6.45e-04): par = (63.8167)
0.06867483  (3.09e-04): par = (63.8458)
0.06867482  (1.48e-04): par = (63.85974)
0.06867481  (7.08e-05): par = (63.86642)
0.06867481  (3.39e-05): par = (63.86962)
0.06867481  (1.62e-05): par = (63.87115)
0.06867481  (7.75e-06): par = (63.87188)
0.08419324  (1.09e-02): par = (66.24142)
0.08417890  (4.73e-03): par = (66.8054)
0.08417621  (2.03e-03): par = (67.05134)
0.08417571  (8.69e-04): par = (67.1572)
0.08417562  (3.71e-04): par = (67.2025)
0.08417561  (1.58e-04): par = (67.22184)
0.08417560  (6.74e-05): par = (67.23009)
0.08417560  (2.87e-05): par = (67.2336)
0.08417560  (1.22e-05): par = (67.2351)
0.08417560  (5.22e-06): par = (67.23574)
0.04219456  (1.98e-02): par = (60.24948)
0.04216898  (1.08e-02): par = (60.92638)
0.04216141  (5.84e-03): par = (61.29767)
0.04215920  (3.14e-03): par = (61.49902)
0.04215856  (1.68e-03): par = (61.60751)
0.04215838  (9.02e-04): par = (61.66577)
0.04215833  (4.83e-04): par = (61.697)
0.04215831  (2.58e-04): par = (61.71371)
0.04215831  (1.38e-04): par = (61.72266)
0.04215830  (7.39e-05): par = (61.72745)
0.04215830  (3.95e-05): par = (61.73001)
0.04215830  (2.11e-05): par = (61.73138)
0.04215830  (1.13e-05): par = (61.73211)
0.04215830  (6.05e-06): par = (61.7325)
0.04729500  (1.53e-02): par = (63.18549)
0.04727818  (8.11e-03): par = (63.76336)
0.04727344  (4.27e-03): par = (64.0723)
0.04727213  (2.24e-03): par = (64.2355)
0.04727176  (1.17e-03): par = (64.32117)
0.04727167  (6.11e-04): par = (64.366)
0.04727164  (3.19e-04): par = (64.3894)
0.04727163  (1.66e-04): par = (64.40162)
0.04727163  (8.67e-05): par = (64.40798)
0.04727163  (4.51e-05): par = (64.4113)
0.04727163  (2.36e-05): par = (64.41303)
0.04727163  (1.23e-05): par = (64.41394)
0.04727163  (6.38e-06): par = (64.41441)
0.05558843  (1.51e-02): par = (67.16139)
0.05556976  (7.28e-03): par = (67.81372)
0.05556541  (3.48e-03): par = (68.13086)
0.05556442  (1.65e-03): par = (68.28277)
0.05556419  (7.83e-04): par = (68.35501)
0.05556414  (3.70e-04): par = (68.38925)
0.05556413  (1.75e-04): par = (68.40545)
0.05556413  (8.27e-05): par = (68.41311)
0.05556413  (3.91e-05): par = (68.41673)
0.05556413  (1.85e-05): par = (68.41844)
0.05556413  (8.74e-06): par = (68.41925)
0.06694045  (1.08e-02): par = (72.57477)
0.06692958  (4.34e-03): par = (73.11744)
0.06692782  (1.73e-03): par = (73.33702)
0.06692754  (6.88e-04): par = (73.4248)
0.06692750  (2.73e-04): par = (73.45972)
0.06692749  (1.08e-04): par = (73.47358)
0.06692749  (4.30e-05): par = (73.47908)
0.06692749  (1.70e-05): par = (73.48126)
0.06692749  (6.76e-06): par = (73.48212)
0.08133989  (6.34e-03): par = (78.79758)
0.08133559  (2.02e-03): par = (79.17211)
0.08133515  (6.43e-04): par = (79.29208)
0.08133510  (2.04e-04): par = (79.33022)
0.08133510  (6.46e-05): par = (79.34231)
0.08133510  (2.05e-05): par = (79.34614)
0.08133510  (6.47e-06): par = (79.34736)
0.09878147  (5.46e-03): par = (85.54564)
0.09877779  (1.37e-03): par = (85.92528)
0.09877756  (3.44e-04): par = (86.02115)
0.09877754  (8.61e-05): par = (86.04519)
0.09877754  (2.16e-05): par = (86.05121)
0.09877754  (5.38e-06): par = (86.05271)
0.04123480  (1.05e-02): par = (390.5248)
0.04123038  (2.50e-04): par = (390.4205)
0.04123037  (5.91e-06): par = (390.423)
0.04624422  (1.16e-02): par = (390.28)
0.04623752  (8.00e-04): par = (390.3974)
0.04623749  (5.48e-05): par = (390.4054)
0.04623749  (3.75e-06): par = (390.406)
0.05669790  (1.41e-02): par = (390.391)
0.05668351  (3.94e-03): par = (390.5488)
0.05668239  (1.10e-03): par = (390.593)
0.05668230  (3.08e-04): par = (390.6053)
0.05668229  (8.59e-05): par = (390.6088)
0.05668229  (2.40e-05): par = (390.6097)
0.05668229  (6.70e-06): par = (390.61)
0.07240720  (2.30e-02): par = (390.7764)
0.07234090  (1.68e-02): par = (391.0942)
0.07230573  (1.21e-02): par = (391.3256)
0.07228772  (8.53e-03): par = (391.491)
0.07227876  (5.96e-03): par = (391.6075)
0.07227442  (4.12e-03): par = (391.6887)
0.07227234  (2.83e-03): par = (391.7447)
0.07227137  (1.93e-03): par = (391.7831)
0.07227091  (1.32e-03): par = (391.8093)
0.07227070  (8.96e-04): par = (391.8271)
0.07227061  (6.08e-04): par = (391.8393)
0.07227056  (4.12e-04): par = (391.8475)
0.07227054  (2.79e-04): par = (391.8531)
0.07227053  (1.89e-04): par = (391.8568)
0.07227053  (1.28e-04): par = (391.8594)
0.07227052  (8.65e-05): par = (391.8611)
0.07227052  (5.85e-05): par = (391.8623)
0.07227052  (3.96e-05): par = (391.8631)
0.07227052  (2.68e-05): par = (391.8636)
0.07227052  (1.81e-05): par = (391.864)
0.07227052  (1.22e-05): par = (391.8642)
0.07227052  (8.27e-06): par = (391.8644)
0.08982678  (1.14e-01): par = (392.8918)
0.08743614  (8.92e-02): par = (394.7994)
0.08638093  (3.49e-02): par = (396.0099)
0.08624132  (1.02e-02): par = (396.4223)
0.08622981  (2.78e-03): par = (396.5376)
0.08622896  (7.44e-04): par = (396.5686)
0.08622890  (1.98e-04): par = (396.5769)
0.08622890  (5.28e-05): par = (396.5791)
0.08622890  (1.41e-05): par = (396.5797)
0.08622890  (3.74e-06): par = (396.5798)
0.04113533  (1.52e-02): par = (390.5248)
0.04112658  (1.22e-03): par = (390.6892)
0.04112652  (9.86e-05): par = (390.676)
0.04112652  (7.96e-06): par = (390.677)
0.04488642  (5.79e-03): par = (390.9047)
0.04488508  (6.28e-04): par = (390.9743)
0.04488507  (6.83e-05): par = (390.9668)
0.04488507  (7.43e-06): par = (390.9676)
0.05175819  (4.94e-03): par = (391.2742)
0.05175710  (6.63e-04): par = (391.3421)
0.05175708  (8.91e-05): par = (391.333)
0.05175708  (1.20e-05): par = (391.3342)
0.05175708  (1.61e-06): par = (391.3341)
0.06170759  (4.64e-03): par = (391.721)
0.06170647  (7.35e-04): par = (391.7952)
0.06170644  (1.17e-04): par = (391.7834)
0.06170644  (1.85e-05): par = (391.7853)
0.06170644  (2.93e-06): par = (391.785)
0.07461996  (5.18e-03): par = (392.2626)
0.07461832  (9.35e-04): par = (392.3595)
0.07461826  (1.69e-04): par = (392.342)
0.07461826  (3.05e-05): par = (392.3452)
0.07461826  (5.50e-06): par = (392.3446)
0.09038888  (5.84e-03): par = (392.9428)
0.09038641  (1.16e-03): par = (393.0708)
0.09038631  (2.33e-04): par = (393.0453)
0.09038631  (4.65e-05): par = (393.0504)
0.09038631  (9.28e-06): par = (393.0493)
Code 
CI_pheno_model
            2.5%    97.5%
center 383.05912 398.4712
width   31.95652 101.2457
Code 
# bootstrapping the output
#bootstrap_CI_pheno_model <- nlsBoot(pheno_model,niter=999)
#summary(bootstrap_CI_pheno_model)
#hist(bootstrap_CI_pheno_model$coefboot[,1], breaks=50)
#hist(bootstrap_CI_pheno_model$coefboot[,2], breaks=50)

#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=pheno_mean,ymin=pheno_mean-pheno_se,ymax=pheno_mean+pheno_se),color="blue",shape=15) +
  geom_smooth(data=sub_summary, aes(x=dist,y=pheno_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="blue") +
  ylim(c(-0.01,1.3)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(face="bold",size=12), axis.text=element_text(size=10,color="black"))+
  geom_rect(aes(xmin = CI_pheno_model[1,1], xmax = CI_pheno_model[1,2], ymin = 1.1, ymax = 1.2), fill = NA, color = "blue") + geom_segment(aes(x = summary_pheno$parameters[1,1], y = 1.1, yend = 1.2), color = "blue", lwd = 1)

0.5 plot them all over each other

Code 
#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=pheno_mean,ymin=pheno_mean-pheno_se,ymax=pheno_mean+pheno_se),color="blue",shape=15) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=ancestry_mean,ymin=ancestry_mean-ancestry_se,ymax=ancestry_mean+ancestry_se),color="gray",shape=16) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=mtDNA_mean,ymin=mtDNA_mean-mtDNA_se,ymax=mtDNA_mean+mtDNA_se),color="black",shape=17) +
  geom_smooth(data=sub_summary, aes(x=dist,y=pheno_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="blue") +
    geom_smooth(data=sub_summary, aes(x=dist,y=ancestry_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="gray") +
  geom_smooth(data=sub_summary, aes(x=dist,y=mtDNA_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="black") +
  ylim(c(-0.01,1.5)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(size=12), axis.text=element_text(size=10,color="black"))+
  geom_rect(aes(xmin = CI_pheno_model[1,1], xmax = CI_pheno_model[1,2], ymin = 1.1, ymax = 1.2), fill = NA, color = "blue") + geom_segment(aes(x = summary_pheno$parameters[1,1], y = 1.1, yend = 1.2), color = "blue", lwd = 1) +
  geom_rect(aes(xmin = CI_ancestry_model[1,1], xmax = CI_ancestry_model[1,2], ymin = 1.22, ymax = 1.32), fill = NA, color = "gray") + geom_segment(aes(x = summary_ancestry$parameters[1,1], y = 1.22, yend = 1.32), color = "gray", lwd = 1) +
  geom_rect(aes(xmin = CI_mtDNA_model[1,1], xmax = CI_mtDNA_model[1,2], ymin = 1.34, ymax = 1.44), fill = NA, color = "black") + geom_segment(aes(x = summarymtdna$parameters[1,1], y = 1.34, yend = 1.44), color = "black", lwd = 1)

Code 
#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=pheno_mean,ymin=pheno_mean-pheno_se,ymax=pheno_mean+pheno_se),color="blue",shape=15) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=ancestry_mean,ymin=ancestry_mean-ancestry_se,ymax=ancestry_mean+ancestry_se),color="gray",shape=16) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=mtDNA_mean,ymin=mtDNA_mean-mtDNA_se,ymax=mtDNA_mean+mtDNA_se),color="black",shape=17) +
  geom_smooth(data=sub_summary, aes(x=dist,y=pheno_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="blue") +
    geom_smooth(data=sub_summary, aes(x=dist,y=ancestry_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="gray") +
  geom_smooth(data=sub_summary, aes(x=dist,y=mtDNA_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="black") +
  ylim(c(-0.01,1.5)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(size=12), axis.text=element_text(size=10,color="black"))+
  geom_segment(aes(x = CI_pheno_model[1,1], xend= CI_pheno_model[1,2], y = 1.1), color = "blue", lwd = 1)+
  geom_point(aes(x = summary_pheno$parameters[1,1], y = 1.1), color = "blue", shape=15, cex=2) +
  geom_segment(aes(x = CI_ancestry_model[1,1], xend= CI_ancestry_model[1,2], y = 1.2), color = "gray", lwd = 1)+
  geom_point(aes(x = summary_ancestry$parameters[1,1], y = 1.2), color = "gray", shape=16, cex=2) +
  geom_segment(aes(x = CI_mtDNA_model[1,1], xend= CI_mtDNA_model[1,2], y = 1.3), color = "black", lwd = 1)+
  geom_point(aes(x=summarymtdna$parameters[1,1], y=1.298), colour="black", shape=17, cex=2)

Code 
pp<-#plot
ggplot() +
  geom_pointrange(data=sub_summary,aes(x=dist,y=pheno_mean,ymin=pheno_mean-pheno_se,ymax=pheno_mean+pheno_se),color="blue",shape=15) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=ancestry_mean,ymin=ancestry_mean-ancestry_se,ymax=ancestry_mean+ancestry_se),color="gray",shape=16) +
    geom_pointrange(data=sub_summary,aes(x=dist,y=mtDNA_mean,ymin=mtDNA_mean-mtDNA_se,ymax=mtDNA_mean+mtDNA_se),color="black",shape=17) +
  geom_smooth(data=sub_summary, aes(x=dist,y=pheno_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="blue") +
    geom_smooth(data=sub_summary, aes(x=dist,y=ancestry_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="gray") +
  geom_smooth(data=sub_summary, aes(x=dist,y=mtDNA_mean), method="nls", formula=y~1/(1+exp((4*(x-c))/w)), se=FALSE, method.args=list(start=list(c=max(sub_summary$dist)/2,w=max(sub_summary$dist)/2)), color="black") +
  ylim(c(-0.01,1.5)) +
  xlim(c(0,max(sub_summary$dist))) +
  xlab("Distance (km)") +
  ylab("Hybrid index") +
  theme_classic() +
  theme(axis.title=element_text(size=12), axis.text=element_text(size=10,color="black"))+
  geom_segment(aes(x = CI_pheno_model[1,1], xend= CI_pheno_model[1,2], y = 1.1), color = "blue", lwd = 1)+
  geom_point(aes(x = summary_pheno$parameters[1,1], y = 1.1), color = "blue", shape=15, cex=2) +
  geom_segment(aes(x = CI_ancestry_model[1,1], xend= CI_ancestry_model[1,2], y = 1.2), color = "gray", lwd = 1)+
  geom_point(aes(x = summary_ancestry$parameters[1,1], y = 1.2), color = "gray", shape=16, cex=2) +
  geom_segment(aes(x = CI_mtDNA_model[1,1], xend= CI_mtDNA_model[1,2], y = 1.3), color = "black", lwd = 1)+
  geom_point(aes(x=summarymtdna$parameters[1,1], y=1.298), colour="black", shape=17, cex=2)

#ggsave("~/Desktop/grosbeak.rad/cline.plot.pdf",pp,width=6,height=4,units="in")