Q & K Box 15.4

Feinberg et al. (2014) examined morphological, genetic, and acoustic (call) criteria in four species of leopard frogs (Rana sphenocephala, R. pipiens, R. palustris, and a new species named R. kauffeldi) and an acoustically similar congener R. sylvatica (acoustic criteria only). For 283 museum specimens across the first four species, they measured size (snout-vent length) and 12 other morphological characteristics: head length, head width, eye diameter, tympanum diameter, foot length, eye to naris distance, naris to snout distance, thigh length, internarial distance, interorbital distance, shank length, and dorsal snout angle. Foot length was not recorded for 19 specimens, so these were excluded from the analysis. They also recorded seven call characteristics (call length, call rate, call rise time, call duty cycle, pulse number, pulse rate, and dominant frequency) from 45 frogs in the field across the five species. Call rate and call length were both adjusted based on regressions against temperature to a standard 14°C. We will not cover the genetic analysis here, and we will use LDA to examine the different species in multivariate space and see if the proposed new species was distinguishable morphologically or acoustically.

Rana sphenocephala. Bob Warrick, , via Wikimedia Commons

Rana pipiens. Ryan Hodnett, , via Wikimedia Commons

Rana kauffeldi. Brian R. Curry, , via Wikimedia Commons

The paper is here and the data sets for this example are morphological feinmorph.csv and we also use the acoustic variables feinacoust.csv

Feinberg, J. A., Newman, C. E., Watkins-Colwell, G. J., Schlesinger, M. D., Zarate, B., Curry, B. R., Shaffer, H. B. & Burger, J. (2014). Cryptic diversity in metropolis: confirmation of a new leopard frog species (Anura: Ranidae) from New York City and surrounding Atlantic coast regions. PLoS One, 9, e108213.

Preliminaries

Load required packages

# packages: vegan, car
source("../R/libraries.R")
library(vegan)
library(DAAG)

Read in data files

feinmorph <- read_csv("../data/feinmorph.csv")

New names:Rows: 264 Columns: 24── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): spp
dbl (13): svl, hw, hl, td, ew, tl, sl, fl, end, nsd, iod, ind, dsa
lgl (10): ...15, ...16, ...17, ...18, ...19, ...20, ...21, ...22, ...23, ...24
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

feinacoust<-read_csv("../data/feinacoust.csv")

Rows: 45 Columns: 8── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): spp
dbl (7): frogid, crt, cdc, pn, df, tcl, tcr
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Check linearity and distributions

scatterplotMatrix(~hw+hl+td+ew+tl+sl+fl+end+nsd+iod+ind+dsa,data=feinmorph,diagonal=list(method='boxplot'))

Standardize variables

feinmorph$hw <- scale(feinmorph$hw)
feinmorph$hl <- scale(feinmorph$hl)
feinmorph$td <- scale(feinmorph$td)
feinmorph$ew <- scale(feinmorph$ew)
feinmorph$tl <- scale(feinmorph$tl)
feinmorph$sl <- scale(feinmorph$sl)
feinmorph$fl <- scale(feinmorph$fl)
feinmorph$end <- scale(feinmorph$end)
feinmorph$nsd <- scale(feinmorph$nsd)
feinmorph$iod <- scale(feinmorph$iod)
feinmorph$ind <- scale(feinmorph$ind)
feinmorph$dsa <- scale(feinmorph$dsa)

Check homog variances and covariances

feinmorph.dist <- dist(feinmorph[,-(1:2)])
feinmorph.disp <- betadisper(feinmorph.dist,feinmorph$spp)
permutest(feinmorph.disp)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
           Df Sum Sq Mean Sq    F N.Perm Pr(>F)
Groups      3     11    3.70 0.86    999   0.47
Residuals 260   1116    4.29                   

anova(feinmorph.disp)

Analysis of Variance Table

Response: Distances
           Df Sum Sq Mean Sq F value Pr(>F)
Groups      3     11    3.70    0.86   0.46
Residuals 260   1116    4.29               

LDA with jackknife CV

feinmorphjac.lda <- lda(spp~hw+hl+td+ew+tl+sl+fl+end+nsd+iod+ind+dsa, data=feinmorph, CV=TRUE)
print(feinmorphjac.lda)

$class
  [1] rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk  
 [39] rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rpa  rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk  
 [77] rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rpa  rk   rk   rk   rk   rk   rk   rpa  rk   rk   rk   rk  
[115] rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rpa  rk   rk   rpa  rk   rk   rk   rsph rpi  rsph rk   rk   rk   rsph rpi  rk   rk  
[153] rk   rk   rk   rk   rk   rk   rk   rpi  rpi  rk   rk   rpi  rpi  rpa  rpi  rpa  rpa  rk   rpa  rk   rk   rk   rk   rpa  rk   rk   rpa  rk   rpa  rk   rpa  rpa  rpa  rk   rk   rpi  rpi  rpi 
[191] rk   rpa  rpi  rsph rpi  rpi  rpi  rpi  rk   rpi  rpi  rk   rpa  rpi  rpi  rpa  rpi  rpi  rk   rpi  rpi  rpi  rpi  rpi  rpi  rpi  rpi  rpi  rsph rsph rsph rsph rsph rsph rsph rk   rsph rsph
[229] rsph rsph rsph rsph rsph rk   rk   rsph rk   rk   rk   rk   rsph rk   rk   rsph rk   rk   rk   rsph rsph rsph rsph rsph rk   rsph rsph rsph rk   rsph rk   rsph rsph rk   rsph rsph
Levels: rk rpa rpi rsph

$posterior
          rk      rpa      rpi     rsph
1   9.68e-01 2.58e-02 1.62e-05 6.22e-03
2   9.49e-01 3.38e-02 1.21e-05 1.67e-02
3   8.16e-01 1.84e-01 2.32e-06 4.99e-04
4   9.73e-01 1.78e-02 1.89e-03 6.87e-03
5   9.84e-01 1.40e-02 1.01e-03 6.10e-04
6   9.22e-01 2.44e-02 8.90e-04 5.32e-02
7   9.46e-01 9.19e-03 2.01e-05 4.46e-02
8   9.34e-01 6.26e-02 9.87e-05 3.41e-03
9   9.42e-01 5.66e-02 2.14e-04 1.56e-03
10  9.56e-01 4.14e-02 1.75e-03 6.61e-04
11  9.80e-01 1.29e-02 3.50e-03 3.59e-03
12  9.06e-01 7.05e-02 1.37e-04 2.30e-02
13  5.53e-01 2.97e-01 1.06e-01 4.39e-02
14  9.58e-01 6.49e-04 3.36e-05 4.17e-02
15  9.21e-01 4.15e-02 2.99e-04 3.71e-02
16  8.92e-01 2.83e-02 1.36e-05 8.00e-02
17  7.29e-01 2.69e-02 2.73e-03 2.41e-01
18  9.30e-01 2.38e-02 2.27e-04 4.59e-02
19  8.36e-01 8.48e-02 1.49e-03 7.80e-02
20  6.29e-01 5.24e-03 1.76e-02 3.48e-01
21  9.77e-01 1.93e-02 6.41e-06 3.20e-03
22  6.71e-01 8.37e-04 5.16e-06 3.29e-01
23  9.97e-01 1.52e-04 3.55e-08 3.18e-03
24  9.66e-01 3.35e-02 2.43e-06 8.66e-04
25  9.60e-01 5.51e-04 8.73e-03 3.02e-02
26  9.80e-01 1.54e-03 1.10e-05 1.82e-02
27  8.62e-01 7.20e-02 1.71e-04 6.60e-02
28  8.92e-01 9.69e-02 3.50e-04 1.11e-02
29  9.23e-01 4.38e-02 1.68e-02 1.62e-02
30  9.95e-01 1.26e-03 1.78e-06 3.41e-03
31  7.06e-01 1.05e-02 2.79e-01 4.26e-03
32  9.46e-01 3.55e-02 1.87e-04 1.88e-02
33  9.82e-01 6.76e-03 3.51e-05 1.11e-02
34  9.74e-01 1.81e-02 8.62e-05 7.98e-03
35  9.89e-01 7.29e-03 1.34e-05 4.17e-03
36  9.51e-01 3.82e-02 3.00e-03 7.47e-03
37  9.29e-01 7.03e-02 1.68e-04 3.36e-04
38  9.38e-01 5.85e-02 2.51e-03 7.47e-04
39  8.64e-01 9.36e-02 7.06e-05 4.21e-02
40  6.69e-01 3.18e-01 1.38e-03 1.16e-02
41  7.84e-01 2.15e-01 7.15e-05 5.05e-04
42  8.08e-01 1.89e-01 6.03e-04 2.18e-03
43  7.60e-01 1.83e-01 4.13e-04 5.67e-02
44  9.59e-01 3.19e-03 3.25e-06 3.81e-02
45  9.43e-01 5.38e-02 2.83e-05 3.30e-03
46  7.67e-01 1.38e-01 6.87e-02 2.67e-02
47  9.93e-01 5.05e-03 1.87e-05 1.87e-03
48  9.87e-01 8.64e-03 1.51e-06 4.16e-03
49  1.53e-01 8.44e-01 2.72e-04 1.98e-03
50  9.30e-01 6.58e-02 4.42e-04 3.36e-03
51  9.99e-01 5.52e-05 4.16e-07 6.55e-04
52  9.95e-01 1.97e-03 1.58e-04 2.50e-03
53  8.91e-01 6.88e-02 2.17e-04 3.96e-02
54  9.61e-01 2.35e-02 4.38e-03 1.08e-02
55  9.87e-01 2.42e-03 1.49e-06 1.08e-02
56  9.97e-01 2.36e-03 1.06e-05 2.04e-04
57  9.93e-01 5.83e-03 1.13e-04 1.22e-03
58  9.98e-01 1.91e-03 5.00e-06 1.19e-04
59  8.70e-01 8.03e-02 4.92e-02 7.38e-04
60  9.20e-01 6.79e-02 8.35e-03 3.52e-03
61  9.85e-01 9.21e-03 1.17e-04 5.79e-03
62  8.89e-01 6.41e-02 7.62e-03 3.88e-02
63  9.77e-01 2.29e-02 3.22e-05 5.27e-04
64  9.70e-01 1.98e-02 1.84e-03 8.45e-03
65  9.98e-01 1.62e-03 9.53e-06 6.58e-05
66  9.98e-01 1.37e-03 5.74e-05 3.59e-04
67  9.74e-01 2.08e-02 6.30e-04 4.64e-03
68  9.79e-01 1.94e-02 8.79e-04 3.07e-04
69  8.66e-01 1.31e-01 9.14e-05 2.64e-03
70  9.70e-01 2.93e-02 4.70e-06 9.60e-04
71  9.75e-01 2.20e-02 2.53e-05 3.38e-03
72  9.37e-01 6.18e-02 3.48e-05 1.30e-03
73  8.93e-01 8.77e-02 6.39e-03 1.32e-02
74  7.73e-01 1.49e-01 5.47e-03 7.31e-02
75  9.92e-01 2.68e-03 2.61e-05 5.06e-03
76  6.20e-01 3.10e-02 2.10e-05 3.49e-01
77  8.71e-01 5.80e-02 1.73e-04 7.08e-02
78  9.94e-01 9.59e-04 1.97e-07 4.78e-03
79  6.26e-01 1.73e-01 5.22e-03 1.96e-01
80  8.74e-01 1.74e-02 1.00e-03 1.08e-01
81  9.46e-01 9.74e-03 1.79e-03 4.23e-02
82  8.18e-01 3.69e-02 4.20e-03 1.41e-01
83  9.90e-01 8.32e-03 1.91e-05 1.62e-03
84  9.87e-01 6.55e-03 2.09e-04 5.82e-03
85  9.00e-01 4.01e-02 8.35e-06 5.95e-02
86  9.73e-01 1.67e-02 7.23e-04 9.25e-03
87  9.93e-01 4.09e-03 3.67e-05 3.13e-03
88  9.59e-01 2.29e-02 1.93e-03 1.64e-02
89  8.35e-01 1.06e-01 2.94e-02 2.99e-02
90  9.85e-01 4.00e-03 4.55e-06 1.11e-02
91  9.93e-01 6.94e-03 3.31e-05 1.24e-04
92  9.97e-01 1.47e-03 1.12e-07 1.65e-03
93  9.96e-01 3.19e-03 2.70e-04 4.91e-04
94  8.75e-01 2.10e-02 7.90e-03 9.65e-02
95  8.07e-01 1.34e-01 3.37e-05 5.90e-02
96  9.94e-01 4.44e-03 6.76e-05 1.17e-03
97  9.89e-01 6.96e-04 5.76e-07 1.01e-02
98  9.36e-01 9.17e-03 7.15e-06 5.45e-02
99  9.88e-01 5.26e-03 3.11e-03 3.33e-03
100 9.80e-01 1.85e-02 3.98e-06 1.20e-03
101 6.38e-01 3.14e-01 4.99e-07 4.74e-02
102 8.41e-01 1.38e-01 2.06e-04 2.06e-02
103 7.95e-02 6.48e-01 1.02e-03 2.71e-01
104 6.86e-01 2.23e-01 2.57e-03 8.84e-02
105 5.70e-01 8.41e-03 4.94e-05 4.22e-01
106 7.12e-01 2.55e-02 1.35e-05 2.62e-01
107 7.92e-01 1.70e-02 9.28e-05 1.91e-01
108 8.78e-01 2.21e-02 3.07e-04 9.96e-02
109 8.06e-01 1.07e-02 1.64e-05 1.83e-01
110 1.15e-01 8.83e-01 4.29e-04 1.42e-03
111 9.98e-01 2.41e-03 1.18e-06 4.46e-05
112 8.44e-01 1.34e-01 7.00e-04 2.20e-02
113 1.00e+00 1.41e-04 2.50e-06 2.15e-04
114 9.70e-01 1.10e-02 4.49e-05 1.91e-02
115 9.65e-01 7.57e-03 1.28e-04 2.72e-02
116 9.88e-01 5.41e-03 2.51e-05 7.03e-03
117 9.38e-01 5.24e-02 9.06e-04 8.51e-03
118 5.95e-01 4.03e-01 1.15e-03 1.60e-03
119 8.66e-01 1.13e-01 1.20e-03 1.95e-02
120 6.89e-01 2.79e-01 2.43e-03 3.02e-02
121 9.55e-01 2.83e-02 2.72e-06 1.63e-02
122 8.30e-01 1.49e-01 5.82e-04 2.05e-02
123 8.97e-01 9.73e-02 1.19e-04 5.63e-03
124 6.94e-01 3.00e-01 8.28e-05 6.18e-03
125 7.71e-01 2.12e-01 3.98e-04 1.67e-02
126 7.15e-01 9.48e-02 4.79e-03 1.85e-01
127 9.38e-01 5.56e-02 5.80e-04 6.26e-03
128 8.17e-01 1.38e-01 1.86e-02 2.59e-02
129 5.86e-01 5.61e-02 8.95e-04 3.57e-01
130 8.48e-01 9.41e-02 1.55e-04 5.74e-02
131 6.65e-01 4.98e-02 2.75e-05 2.85e-01
132 7.64e-01 1.98e-01 2.46e-04 3.85e-02
133 9.31e-01 3.82e-02 2.72e-04 3.06e-02
134 9.30e-01 2.01e-02 1.09e-05 5.00e-02
135 7.35e-01 2.63e-01 1.14e-04 2.04e-03
136 1.96e-04 9.99e-01 5.05e-04 1.45e-05
137 7.35e-01 8.92e-03 9.57e-04 2.55e-01
138 9.84e-01 3.98e-03 9.60e-05 1.19e-02
139 2.98e-02 9.69e-01 5.90e-05 1.23e-03
140 9.19e-01 6.89e-02 6.82e-06 1.23e-02
141 6.77e-01 3.08e-01 1.13e-02 3.21e-03
142 9.56e-01 3.24e-02 1.35e-03 1.07e-02
143 1.86e-01 1.39e-02 6.66e-02 7.33e-01
144 3.09e-01 4.55e-02 4.60e-01 1.85e-01
145 2.53e-01 2.94e-02 8.50e-04 7.17e-01
146 6.93e-01 1.63e-02 4.01e-04 2.90e-01
147 4.89e-01 1.11e-01 1.68e-03 3.99e-01
148 4.95e-01 9.53e-02 8.26e-03 4.02e-01
149 1.13e-01 1.21e-02 2.24e-04 8.75e-01
150 1.03e-01 1.70e-01 7.27e-01 2.40e-05
151 9.59e-01 3.84e-02 4.72e-04 1.67e-03
152 8.22e-01 5.16e-02 9.27e-03 1.18e-01
153 8.97e-01 7.44e-02 4.61e-03 2.38e-02
154 9.96e-01 2.23e-03 3.18e-04 1.48e-03
155 8.79e-01 1.24e-02 4.77e-03 1.04e-01
156 9.82e-01 1.41e-02 5.75e-06 3.43e-03
157 8.99e-01 9.74e-02 1.77e-04 3.04e-03
158 7.03e-01 2.68e-01 1.24e-02 1.65e-02
159 6.33e-01 3.58e-01 1.58e-05 8.61e-03
160 4.48e-04 8.39e-04 9.95e-01 4.20e-03
161 3.40e-02 5.60e-02 9.03e-01 7.36e-03
162 4.14e-01 3.60e-01 3.13e-05 2.26e-01
163 8.14e-01 1.32e-01 8.49e-04 5.31e-02
164 6.72e-02 1.26e-01 7.49e-01 5.74e-02
165 9.14e-04 9.71e-03 9.89e-01 5.85e-06
166 4.21e-01 4.74e-01 9.83e-02 7.20e-03
167 3.54e-03 1.99e-01 7.90e-01 7.32e-03
168 1.37e-01 7.13e-01 1.37e-01 1.20e-02
169 2.99e-01 6.15e-01 1.00e-03 8.53e-02
170 7.90e-01 1.61e-01 8.07e-04 4.75e-02
171 1.70e-01 8.05e-01 4.53e-04 2.42e-02
172 6.46e-01 3.49e-01 1.93e-03 2.16e-03
173 9.53e-01 4.53e-02 1.76e-04 1.28e-03
174 8.48e-01 1.14e-01 9.29e-05 3.82e-02
175 8.50e-01 8.88e-02 5.85e-05 6.09e-02
176 4.06e-01 5.69e-01 5.55e-04 2.47e-02
177 9.57e-01 2.52e-02 4.62e-05 1.76e-02
178 8.48e-01 7.03e-02 7.63e-03 7.38e-02
179 3.81e-01 5.78e-01 1.13e-03 3.99e-02
180 5.94e-01 2.22e-01 1.28e-02 1.71e-01
181 1.84e-01 8.04e-01 6.84e-05 1.21e-02
182 6.45e-01 3.27e-01 6.64e-05 2.70e-02
183 1.37e-01 8.61e-01 1.08e-03 1.42e-03
184 2.46e-02 8.54e-01 1.08e-01 1.32e-02
185 2.87e-01 6.84e-01 2.43e-02 4.85e-03
186 8.50e-01 1.47e-01 7.18e-04 2.07e-03
187 9.51e-01 4.24e-02 2.51e-04 6.83e-03
188 7.50e-03 5.79e-03 7.17e-01 2.70e-01
189 1.44e-02 6.00e-02 9.11e-01 1.45e-02
190 3.00e-04 2.22e-04 9.99e-01 4.37e-05
191 8.26e-01 1.58e-01 2.77e-03 1.32e-02
192 2.11e-02 6.35e-01 2.60e-01 8.48e-02
193 1.23e-02 2.84e-03 9.76e-01 9.21e-03
194 1.09e-01 4.01e-03 1.29e-02 8.75e-01
195 1.90e-02 1.74e-03 8.07e-01 1.73e-01
196 5.43e-06 3.59e-05 9.99e-01 1.27e-03
197 4.10e-04 5.94e-04 9.99e-01 1.34e-04
198 5.08e-03 2.94e-02 9.62e-01 3.75e-03
199 3.25e-01 2.98e-01 2.33e-01 1.44e-01
200 1.71e-01 1.88e-01 5.77e-01 6.42e-02
201 3.16e-04 1.32e-03 9.98e-01 2.11e-05
202 4.86e-01 7.78e-02 1.57e-01 2.80e-01
203 7.02e-02 6.00e-01 3.24e-01 5.90e-03
204 2.28e-05 5.41e-04 9.99e-01 3.00e-05
205 1.00e-04 6.85e-05 1.00e+00 8.54e-05
206 2.79e-01 5.55e-01 1.63e-01 3.18e-03
207 3.54e-03 9.17e-03 9.86e-01 1.39e-03
208 2.50e-04 3.33e-03 9.96e-01 9.86e-06
209 8.05e-01 9.97e-02 7.09e-02 2.46e-02
210 3.80e-04 3.37e-04 9.99e-01 9.24e-05
211 7.09e-04 3.36e-03 9.96e-01 1.83e-04
212 3.01e-05 3.72e-03 9.96e-01 4.45e-05
213 3.72e-03 6.91e-03 9.84e-01 5.12e-03
214 2.57e-02 1.50e-03 9.72e-01 5.72e-04
215 2.54e-05 8.98e-04 9.98e-01 9.14e-04
216 4.99e-02 1.78e-01 7.50e-01 2.24e-02
217 1.40e-02 3.96e-02 9.42e-01 4.04e-03
218 1.40e-04 1.21e-03 9.98e-01 4.10e-04
219 5.03e-03 1.68e-02 4.71e-05 9.78e-01
220 1.28e-04 5.97e-05 2.47e-01 7.53e-01
221 2.01e-03 9.86e-04 2.41e-01 7.56e-01
222 2.28e-07 3.03e-07 4.37e-08 1.00e+00
223 5.64e-02 1.48e-01 1.11e-01 6.84e-01
224 1.26e-05 2.59e-06 1.73e-07 1.00e+00
225 3.08e-02 1.77e-03 1.17e-02 9.56e-01
226 4.92e-01 9.67e-02 3.73e-02 3.74e-01
227 7.56e-02 3.65e-03 5.69e-04 9.20e-01
228 3.15e-02 5.33e-03 2.93e-03 9.60e-01
229 1.74e-02 5.96e-02 4.10e-02 8.82e-01
230 1.18e-01 4.84e-02 2.16e-02 8.12e-01
231 1.21e-02 2.21e-02 3.54e-02 9.30e-01
232 2.58e-02 5.61e-03 2.73e-04 9.68e-01
233 2.78e-02 5.41e-03 3.59e-03 9.63e-01
234 6.78e-01 5.71e-02 2.58e-03 2.62e-01
235 6.82e-01 3.65e-02 3.08e-04 2.82e-01
236 1.17e-01 2.70e-03 3.29e-06 8.80e-01
237 5.99e-01 1.80e-01 1.31e-01 9.03e-02
238 7.47e-01 1.64e-01 5.34e-03 8.36e-02
239 5.65e-01 1.50e-01 8.36e-04 2.84e-01
240 9.52e-01 6.16e-03 2.72e-05 4.23e-02
241 4.26e-01 9.31e-02 2.74e-04 4.81e-01
242 9.75e-01 2.06e-02 8.15e-06 4.78e-03
243 9.05e-01 8.86e-02 2.55e-04 6.61e-03
244 6.98e-02 7.15e-02 4.63e-05 8.59e-01
245 9.49e-01 4.98e-02 1.63e-05 7.50e-04
246 8.46e-01 1.37e-01 3.56e-04 1.62e-02
247 5.20e-01 2.78e-01 2.09e-02 1.81e-01
248 1.47e-01 5.76e-03 3.42e-02 8.13e-01
249 3.43e-02 6.67e-03 1.28e-02 9.46e-01
250 4.37e-03 5.51e-04 3.29e-05 9.95e-01
 [ reached getOption("max.print") -- omitted 14 rows ]

$terms
spp ~ hw + hl + td + ew + tl + sl + fl + end + nsd + iod + ind + 
    dsa
attr(,"variables")
list(spp, hw, hl, td, ew, tl, sl, fl, end, nsd, iod, ind, dsa)
attr(,"factors")
    hw hl td ew tl sl fl end nsd iod ind dsa
spp  0  0  0  0  0  0  0   0   0   0   0   0
hw   1  0  0  0  0  0  0   0   0   0   0   0
hl   0  1  0  0  0  0  0   0   0   0   0   0
td   0  0  1  0  0  0  0   0   0   0   0   0
ew   0  0  0  1  0  0  0   0   0   0   0   0
tl   0  0  0  0  1  0  0   0   0   0   0   0
sl   0  0  0  0  0  1  0   0   0   0   0   0
fl   0  0  0  0  0  0  1   0   0   0   0   0
end  0  0  0  0  0  0  0   1   0   0   0   0
nsd  0  0  0  0  0  0  0   0   1   0   0   0
iod  0  0  0  0  0  0  0   0   0   1   0   0
ind  0  0  0  0  0  0  0   0   0   0   1   0
dsa  0  0  0  0  0  0  0   0   0   0   0   1
attr(,"term.labels")
 [1] "hw"  "hl"  "td"  "ew"  "tl"  "sl"  "fl"  "end" "nsd" "iod" "ind" "dsa"
attr(,"order")
 [1] 1 1 1 1 1 1 1 1 1 1 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
attr(,"predvars")
list(spp, hw, hl, td, ew, tl, sl, fl, end, nsd, iod, ind, dsa)
attr(,"dataClasses")
        spp          hw          hl          td          ew          tl          sl          fl         end         nsd         iod         ind         dsa 
"character" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" 

$call
lda(formula = spp ~ hw + hl + td + ew + tl + sl + fl + end + 
    nsd + iod + ind + dsa, data = feinmorph, CV = TRUE)

$xlevels
named list()

table(feinmorph$spp,feinmorphjac.lda$class)

      
        rk rpa rpi rsph
  rk   147   5   2    3
  rpa   15  10   5    0
  rpi    4   3  23    1
  rsph  16   0   0   30

Use Borcard code to get props correct

feinmorphjac.class <- feinmorphjac.lda$class
feinmorphjac.table <- table(feinmorph$spp, feinmorphjac.class)
diag(prop.table(feinmorphjac.table,1))

   rk   rpa   rpi  rsph 
0.936 0.333 0.742 0.652 

Use DAAG confusion to get overall accuracy

confusion(feinmorphjac.lda$class,feinmorph$spp)

Overall accuracy = 0.795 

Confusion matrix 
      Predicted (cv)
Actual    rk   rpa   rpi  rsph
  rk   0.808 0.082 0.022 0.088
  rpa  0.278 0.556 0.167 0.000
  rpi  0.067 0.167 0.767 0.000
  rsph 0.088 0.000 0.029 0.882

LDA

feinmorph.lda <- lda(spp~hw+hl+td+ew+tl+sl+fl+end+nsd+iod+ind+dsa, data=feinmorph)
print(feinmorph.lda)

Call:
lda(spp ~ hw + hl + td + ew + tl + sl + fl + end + nsd + iod + 
    ind + dsa, data = feinmorph)

Prior probabilities of groups:
   rk   rpa   rpi  rsph 
0.595 0.114 0.117 0.174 

Group means:
          hw      hl      td     ew     tl     sl     fl    end    nsd    iod      ind     dsa
rk    0.0837 -0.0326  0.1971 -0.247 -0.184 -0.112 -0.006 -0.244 -0.144  0.280  0.00254  0.2006
rpa  -0.3719 -0.4452 -0.7933 -0.600 -0.211 -0.302 -0.520 -0.221 -0.438 -0.421 -0.19299  0.0297
rpi   0.1812 -0.0542 -0.2849  1.094  0.578  0.590  0.284  0.345  0.923 -0.462  0.52526  0.4524
rsph -0.1652  0.4381  0.0367  0.497  0.375  0.181  0.168  0.744  0.156 -0.368 -0.23680 -1.0089

Coefficients of linear discriminants:
         LD1      LD2     LD3
hw   0.00512  3.21961  3.0955
hl  -1.39156 -3.19366 -3.2937
td  -0.88528 -0.30354  0.6674
ew   0.70656 -0.13599  0.9049
tl   0.16382 -0.14807 -0.7371
sl   2.01583  1.23039 -0.8579
fl  -0.28228 -0.51123  1.2329
end  0.03603 -0.88537 -0.0415
nsd  0.45803  0.14663  0.5644
iod -0.45060 -0.00645  0.2566
ind -0.06066  0.53421 -0.3457
dsa -0.25463 -0.89680 -1.6415

Proportion of trace:
   LD1    LD2    LD3 
0.5959 0.3044 0.0996 

# get predictions from this LDA
feinmorph.pred <- predict(feinmorph.lda)
table(feinmorph$spp,feinmorph.pred$class)

      
        rk rpa rpi rsph
  rk   148   5   1    3
  rpa   14  11   5    0
  rpi    3   2  25    1
  rsph  14   0   0   32

# plot LDA
plot(feinmorph.lda, dimen=2)

Now analyse acoustic data as for morphology

Diagnostics

# check linearity and distributions
scatterplotMatrix(~crt+cdc+pn+df+tcl+tcr,data=feinacoust,diagonal=list(method='boxplot'))

# standardize variables
feinacoust$crt <- scale(feinacoust$crt)
feinacoust$cdc <- scale(feinacoust$cdc)
feinacoust$pn <- scale(feinacoust$pn)
feinacoust$df <- scale(feinacoust$df)
feinacoust$tcl <- scale(feinacoust$tcl)
feinacoust$tcr <- scale(feinacoust$tcr)
# check var-cov homogeneity
feinacoust.dist <- dist(feinacoust[,-(1:2)])
feinacoust.disp <- betadisper(feinacoust.dist,feinacoust$spp)
permutest(feinacoust.disp)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq Mean Sq    F N.Perm Pr(>F)  
Groups     4   5.93   1.482 4.02    999  0.014 *
Residuals 40  14.73   0.368                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

LDA with jackknife CV

feinacoustjac.lda <- lda(spp~crt+cdc+pn+df+tcl+tcr, CV=TRUE, data=feinacoust)
print(feinacoustjac.lda)

$class
 [1] rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rk   rp   rp   rp   rp   rsph rsph rsph rsph rsph rsph rsph rsph rsyl rsyl rsyl rsyl rk   rsyl rsyl rk   rsyl rpal rpal rpal rpal
[39] rpal rpal rpal rpal rpal rpal rpal
Levels: rk rp rpal rsph rsyl

$posterior
          rk       rp     rpal     rsph     rsyl
1   1.00e+00 1.74e-43 6.24e-36 2.12e-15 3.50e-06
2   1.00e+00 2.05e-42 7.24e-36 1.53e-14 1.00e-05
3   1.00e+00 1.35e-42 2.18e-36 1.54e-14 6.15e-05
4   1.00e+00 2.36e-43 2.17e-35 4.67e-15 1.32e-06
5   1.00e+00 1.84e-42 8.32e-36 1.84e-14 6.47e-06
6   1.00e+00 5.65e-44 1.95e-35 3.74e-15 1.70e-06
7   1.00e+00 2.62e-42 2.25e-35 3.64e-14 4.92e-06
8   1.00e+00 2.54e-42 1.85e-35 6.17e-14 5.60e-06
9   1.00e+00 1.17e-45 1.46e-35 6.56e-16 4.54e-06
10  1.00e+00 3.20e-45 5.93e-36 8.16e-16 1.23e-05
11  1.00e+00 2.09e-44 3.36e-35 1.71e-15 1.04e-06
12  1.00e+00 1.44e-39 2.29e-33 9.36e-12 4.47e-07
13  1.00e+00 4.31e-43 4.98e-35 2.44e-14 2.97e-06
14  2.24e-39 1.00e+00 4.48e-33 2.94e-26 5.47e-32
15  6.80e-33 1.00e+00 1.46e-22 1.75e-21 1.45e-25
16  1.30e-33 1.00e+00 2.88e-28 1.91e-19 2.25e-25
17 3.48e-129 1.00e+00 1.85e-55 1.99e-74 1.56e-90
18  4.66e-10 1.05e-30 2.76e-33 1.00e+00 3.53e-09
19  3.56e-13 5.79e-28 2.94e-35 1.00e+00 1.14e-11
20  2.03e-17 4.06e-27 1.60e-36 1.00e+00 7.62e-15
21  3.44e-16 2.45e-27 1.56e-40 1.00e+00 7.16e-15
22  2.70e-17 2.25e-25 2.72e-32 1.00e+00 2.32e-13
23  3.80e-09 1.18e-28 4.05e-29 1.00e+00 1.05e-07
24  1.19e-13 1.37e-29 1.65e-31 1.00e+00 2.15e-10
25  4.97e-18 1.87e-28 1.93e-33 1.00e+00 3.18e-13
26  1.02e-03 3.39e-26 1.08e-31 2.83e-01 7.16e-01
27  5.42e-08 2.31e-32 1.80e-38 1.66e-08 1.00e+00
28  1.33e-07 9.55e-43 2.73e-41 8.13e-14 1.00e+00
29  3.93e-10 4.40e-38 1.19e-40 2.00e-10 1.00e+00
30  9.39e-01 4.39e-31 6.57e-33 2.72e-09 6.06e-02
31  4.90e-09 5.88e-48 4.17e-52 2.58e-21 1.00e+00
32  9.32e-04 1.35e-25 1.47e-35 1.83e-14 9.99e-01
33  9.97e-01 1.55e-32 8.34e-34 3.40e-12 2.86e-03
34  1.91e-10 2.90e-40 1.23e-49 1.31e-27 1.00e+00
35  2.68e-23 6.40e-22 1.00e+00 3.81e-22 8.65e-27
36  3.08e-36 1.82e-23 1.00e+00 7.90e-30 2.85e-37
37  1.50e-44 1.19e-40 1.00e+00 1.60e-40 2.22e-47
38  1.54e-28 3.09e-35 1.00e+00 1.66e-27 1.24e-32
39  5.62e-36 5.26e-38 1.00e+00 1.64e-32 1.29e-40
40  3.81e-36 3.45e-40 1.00e+00 1.81e-32 2.25e-39
41  2.64e-12 1.26e-08 1.00e+00 1.31e-13 5.18e-15
42  1.42e-53 2.30e-53 1.00e+00 2.46e-64 1.90e-62
43  1.34e-34 1.23e-29 1.00e+00 1.76e-34 6.09e-38
44  1.62e-71 6.00e-42 1.00e+00 2.36e-68 9.43e-77
45  5.49e-34 9.81e-33 1.00e+00 3.14e-36 1.37e-37

$terms
spp ~ crt + cdc + pn + df + tcl + tcr
attr(,"variables")
list(spp, crt, cdc, pn, df, tcl, tcr)
attr(,"factors")
    crt cdc pn df tcl tcr
spp   0   0  0  0   0   0
crt   1   0  0  0   0   0
cdc   0   1  0  0   0   0
pn    0   0  1  0   0   0
df    0   0  0  1   0   0
tcl   0   0  0  0   1   0
tcr   0   0  0  0   0   1
attr(,"term.labels")
[1] "crt" "cdc" "pn"  "df"  "tcl" "tcr"
attr(,"order")
[1] 1 1 1 1 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
attr(,"predvars")
list(spp, crt, cdc, pn, df, tcl, tcr)
attr(,"dataClasses")
        spp         crt         cdc          pn          df         tcl         tcr 
"character" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" "nmatrix.1" 

$call
lda(formula = spp ~ crt + cdc + pn + df + tcl + tcr, data = feinacoust, 
    CV = TRUE)

$xlevels
named list()

table(feinacoust$spp,feinacoustjac.lda$class)

      
       rk rp rpal rsph rsyl
  rk   13  0    0    0    0
  rp    0  4    0    0    0
  rpal  0  0   11    0    0
  rsph  0  0    0    8    0
  rsyl  2  0    0    0    7

# use Borcard code to get props correct
feinacoustjac.class <- feinacoustjac.lda$class
feinacoustjac.table <- table(feinacoust$spp, feinacoustjac.class)
diag(prop.table(feinacoustjac.table,1))

   rk    rp  rpal  rsph  rsyl 
1.000 1.000 1.000 1.000 0.778 

# use DAAG confusion to get overall accuracy
confusion(feinacoustjac.lda$class,feinacoust$spp)

Overall accuracy = 0.956 

Confusion matrix 
      Predicted (cv)
Actual    rk rp rpal rsph  rsyl
  rk   0.867  0    0    0 0.133
  rp   0.000  1    0    0 0.000
  rpal 0.000  0    1    0 0.000
  rsph 0.000  0    0    1 0.000
  rsyl 0.000  0    0    0 1.000

LDA

feinacoust.lda <- lda(spp~crt+cdc+pn+df+tcl+tcr, data=feinacoust)
print(feinacoust.lda)

Call:
lda(spp ~ crt + cdc + pn + df + tcl + tcr, data = feinacoust)

Prior probabilities of groups:
    rk     rp   rpal   rsph   rsyl 
0.2889 0.0889 0.2444 0.1778 0.2000 

Group means:
         crt     cdc     pn     df    tcl     tcr
rk   -0.9337 -0.9346 -0.754  0.329 -1.005  0.0629
rp    1.9432 -0.6713  0.703 -0.668  2.117 -0.9856
rpal  0.9377 -0.0902  1.562 -0.240  0.898 -0.8941
rsph -0.0464  1.6022 -0.490 -0.477 -0.282 -0.1782
rsyl -0.6198  0.3344 -0.696  0.538 -0.336  1.5983

Coefficients of linear discriminants:
        LD1    LD2    LD3    LD4
crt  1.1811 -0.261 -0.561  2.004
cdc -0.2784  1.104  2.283 -0.830
pn  -4.9959 -5.147  0.573 -1.389
df   0.0505 -0.525 -0.180 -0.138
tcl -1.1535  5.497 -1.504 -1.691
tcr  0.8200 -1.065 -1.687 -1.463

Proportion of trace:
   LD1    LD2    LD3    LD4 
0.6118 0.2386 0.1063 0.0433 

# get eigenvalues (prop explained)
feinacoust.lda$svd

[1] 17.97 11.22  7.49  4.78

propexp <- feinacoust.lda$svd^2/sum(feinacoust.lda$svd^2)
propexp

[1] 0.6118 0.2386 0.1063 0.0433

feinacoust.pred <- predict(feinacoust.lda)
table(feinacoust$spp,feinacoust.pred$class)

      
       rk rp rpal rsph rsyl
  rk   13  0    0    0    0
  rp    0  4    0    0    0
  rpal  0  0   11    0    0
  rsph  0  0    0    8    0
  rsyl  2  0    0    0    7

plot(feinacoust.lda, dimen=2)

Finally, do MANOVAs

morph <- cbind(feinmorph$hw,feinmorph$hl,feinmorph$td,feinmorph$ew,feinmorph$tl,feinmorph$sl,feinmorph$fl,feinmorph$end,feinmorph$nsd,feinmorph$iod,feinmorph$ind,feinmorph$dsa)
feinmorph.man <- manova(morph~spp, data=feinmorph)
summary(feinmorph.man, test="Pillai")

           Df Pillai approx F num Df den Df Pr(>F)    
spp         3    1.2     13.9     36    753 <2e-16 ***
Residuals 260                                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

acoust <- cbind(feinacoust$crt,feinacoust$cdc,feinacoust$pn,feinacoust$df,feinacoust$tcl,feinacoust$tcr)
feinacoust.man <- manova(acoust~spp, data=feinacoust)
summary(feinacoust.man, test="Pillai")

          Df Pillai approx F num Df den Df Pr(>F)    
spp        4   3.44       39     24    152 <2e-16 ***
Residuals 40                                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Better graphics

Plots used for QK use the ggplot classic theme, with some tweaks. Tweaks are consolidated into theme_QK; use this theme for figures and tweak the theme to avoid repetitive code changes.

source("../R/appearance.R")

lda.data <- cbind(feinmorph, predict(feinmorph.lda)$x)
#Get complex hulls
hull <- 
  lda.data %>%
  group_by(spp) %>% 
  slice(chull(LD1, LD2))
la=c("R.kauffeldi", "R. pipiens", "R. palustris", "R. spheno.", "R. sylvatica")
p1<-ggplot(lda.data, aes(LD1, LD2)) +
  geom_point(aes(shape = spp), show.legend=FALSE, size= ss/2)+
  scale_shape_manual(values=c(0:5))+
  geom_polygon(data=hull, aes(fill = spp,
                   ), color=lc,
               alpha = 0.3,
               show.legend = FALSE)+
  scale_fill_grey(start = 0.4, end = .9)
  
lda2.data <- cbind(feinacoust, predict(feinacoust.lda)$x)
#Get complex hulls
hull2 <- 
  lda2.data %>%
  drop_na() %>%
  group_by(spp) %>% 
  slice(chull(LD1, LD2))
p2<-ggplot(lda2.data, aes(LD1, LD2)) +
  geom_point(aes(shape = spp), size=ss/2)+
  scale_shape_manual(values=c(0:5),
                     name="Species",
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(angle = 0, face = "italic", size=6),
                                    title=NULL))+
  labs(y=NULL)+
  geom_polygon(data=hull2, aes(fill = spp,
                   ), colour = "black",
               alpha = 0.3,
               show.legend = FALSE)+
  scale_fill_grey(start = 0.4, end = .9)
   

p3<-p1+p2 & theme_qk()

p3

Colour fig

p1a<-ggplot(lda.data, aes(LD1, LD2, colour=spp)) +
  geom_point(aes(), show.legend=FALSE, size= ss/2)+
  scale_color_viridis_d()+
  geom_polygon(data=hull, aes(fill = spp,
                   ), 
               alpha = 0.3,
               show.legend = FALSE)+
  scale_fill_viridis_d()
  
p2a<-ggplot(lda2.data, aes(LD1, LD2, colour=spp)) +
  geom_point(size=ss/2)+
  scale_color_viridis_d(
                     name="Species",
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(angle = 0, face = "italic", size=6),
                                    title=NULL))+
  labs(y=NULL)+
  geom_polygon(data=hull2, aes(fill = spp,
                   ),
               alpha = 0.3,
               show.legend = FALSE)+
  scale_fill_viridis_d()
p3a<-p1a+p2a&theme_qk()   

p3a