par.opt: Parameter optimization
In SongEvo: An Individual-Based Model of Bird Song Evolution
Description
Usage
Arguments
Value
See Also
Examples
Description
This function follows par.sens to help users optimize values for imperfectly known parameters for SongEvo. The goals are to maximize accuracy and precision of model prediction.
Usage
1
Arguments
sens.results
The sens.results array from par.sens(), which includes summary.results from SongEvo() for a range of parameter values. summary.results from SongEvo() includes population summary values for each time step (dimension 1) in each iteration (dimension 2) of the model. Population summary values are contained in five additional dimensions: population size for each time step of each iteration (<e2><80><9c>sample.n<e2><80><9d>), the population mean and variance of the song feature studied (<e2><80><9c>trait.pop.mean<e2><80><9d> and <e2><80><9c>trait.pop.variance<e2><80><9d>), with associated lower (<e2><80><9c>lci<e2><80><9d>) and upper (<e2><80><9c>uci<e2><80><9d>) confidence intervals.
ts
The timestep (<e2><80><9c>ts<e2><80><9d>) at which to compare simulated trait values to target trait values (<e2><80><9c>target.data<e2><80><9d>).
target.data
A set of trait values against which to compare simulated values. target.data may be measured (e.g. from a training population) or hypothetical.
par.range
Range of parameter values over which to optimize.
Value
Three measurements of accuracy and one measure of precision. Accuracy is quantified by three different approaches: i) the mean of absolute residuals of the predicted population mean values in relation to observed values (smaller absolute residuals indicate a more accurate model), ii) the difference between the bootstrapped mean of predicted population means and the mean of the observed values, and iii) the proportion of simulated population trait means that fall within confidence intervals of the observed data (a higher proportion indicates greater accuracy). Precision is measured with the residuals of the predicted population variance to the variance of observed values (smaller residuals indicate a more precise model).
See Also
SongEvo, par.sens, mod.val, h.test
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
### See vignette for an example that uses all functions in SongEvo.
#### Specify and call `par.sens()`
# Here we test the sensitivity of the Acquire a Territory submodel to variation
# in territory turnover rates, ranging from 0.8<e2><80><93>1.2 times the published rate
# (40<e2><80><93>60% of territories turned over). The call for the par.sens function has a
# format similar to SongEvo. The user specifies the parameter to test and the
# range of values for that parameter. The function currently allows examination
# of only one parameter at a time and requires at least two iterations.
parm <- "terr.turnover"
par.range = seq(from=0.45, to=0.55, by=0.05)
sens.results <- NULL
data("song.data")
data("glo.parms")
years=2005-1969
iteration=5
timestep=1
n.territories <- glo.parms$n.territories
starting.trait <- subset(song.data, Population=="PRBO" & Year==1969)$Trill.FBW
starting.trait2 <- c(starting.trait, rnorm(n.territories-length(starting.trait),
mean=mean(starting.trait), sd=sd(starting.trait)))
init.inds <- data.frame(id = seq(1:n.territories), age = 2, trait = starting.trait2)
init.inds$x1 <- round(runif(n.territories, min=-122.481858, max=-122.447270), digits=8)
init.inds$y1 <- round(runif(n.territories, min=37.787768, max=37.805645), digits=8)
# Now we call the par.sens function with our specifications.
extra_parms <- list(init.inds = init.inds,
timestep = 1,
n.territories = nrow(init.inds),
learning.method = "integrate",
integrate.dist = 0.1,
lifespan = NA,
terr.turnover = 0.5,
mate.comp = FALSE,
prin = FALSE,
all = TRUE)
global_parms_key <- which(!names(glo.parms) %in% names(extra_parms))
extra_parms[names(glo.parms[global_parms_key])]=glo.parms[global_parms_key]
par.sens1 <- par.sens(parm = parm, par.range = par.range,
iteration = iteration, steps = years, mate.comp = FALSE,
fixed_parms=extra_parms[names(extra_parms)!=parm], all = TRUE)
#### Prepare current song values
target.data <- subset(song.data, Population=="PRBO" & Year==2005)$Trill.FBW
#### Specify and call `par.opt()`
# Users specify the timestep (<e2><80><9c>ts<e2><80><9d>) at which to compare simulated trait values
# to target trait data (<e2><80><9c>target.data<e2><80><9d>) and save the results in an object (called
# `par.opt1` here).
ts <- years
par.opt1 <- par.opt(sens.results=par.sens1$sens.results, ts=ts,
target.data=target.data, par.range=par.range)
# Examine results objects (residuals and target match).
par.opt1$Residuals
par.opt1$Target.match
#### Plot results of `par.opt()`
#### Accuracy
#1. Difference in means.
plot(par.range, par.opt1$Target.match[,1], type="l",
xlab="Parameter range", ylab="Difference in means (Hz)")
#2. Plot proportion contained.
plot(par.range, par.opt1$Prop.contained, type="l",
xlab="Parameter range", ylab="Proportion contained")
#3. Calculate and plot mean and quantiles of residuals of mean trait values.
res.mean.means <- apply(par.opt1$Residuals[, , 1], MARGIN=1,
mean, na.rm=TRUE)
res.mean.quants <- apply (par.opt1$Residuals[, , 1], MARGIN=1,
quantile, probs=c(0.975, 0.025), R=600, na.rm=TRUE)
plot(par.range, res.mean.means, col="orange",
ylim=range(par.opt1$Residuals[,,1], na.rm=TRUE),
type="b",
xlab="Parameter value (territory turnover rate)",
ylab="Residual of trait mean (trill bandwidth, Hz)")
points(par.range, res.mean.quants[1,], col="orange")
points(par.range, res.mean.quants[2,], col="orange")
lines(par.range, res.mean.quants[1,], col="orange", lty=2)
lines(par.range, res.mean.quants[2,], col="orange", lty=2)
#### Precision
#Calculate and plot mean and quantiles of residuals of variance of trait values
res.var.mean <- apply(par.opt1$Residuals[, , 2], MARGIN=1,
mean, na.rm=TRUE)
res.var.quants <- apply (par.opt1$Residuals[, , 2], MARGIN=1,
quantile, probs=c(0.975, 0.025), R=600, na.rm=TRUE)
plot(par.range, res.var.mean, col="purple",
ylim=range(par.opt1$Residuals[,,2], na.rm=TRUE),
type="b",
xlab="Parameter value (territory turnover rate)",
ylab="Residual of trait variance (trill bandwidth, Hz)")
points(par.range, res.var.quants[1,], col="purple")
points(par.range, res.var.quants[2,], col="purple")
lines(par.range, res.var.quants[1,], col="purple", lty=2)
lines(par.range, res.var.quants[2,], col="purple", lty=2)
#### Visual inspection of accuracy and precision: plot trait values for range of parameters
par(mfcol=c(3,2),
mar=c(4.1, 4.1, 1, 1),
cex=1.2)
for(i in 1:length(par.range)){
plot(par.sens1$sens.results[ , , "trait.pop.mean", ],
xlab="Year", ylab="Bandwidth (Hz)",
xaxt="n", type="n",
xlim=c(-0.5, years), ylim=range(par.sens1$sens.results[ , , "trait.pop.mean", ], na.rm=TRUE))
for(p in 1:iteration){
lines(par.sens1$sens.results[p, , "trait.pop.mean", i], col="light gray")
}
freq.mean <- apply(par.sens1$sens.results[, , "trait.pop.mean", i], 2, mean, na.rm=TRUE)
lines(freq.mean, col="blue")
axis(side=1, at=seq(0, 35, by=5), labels=seq(1970, 2005, by=5))#, tcl=-0.25, mgp=c(2,0.5,0))
#Plot 95% quantiles
quant.means <- apply (par.sens1$sens.results[, , "trait.pop.mean", i], MARGIN=2,
quantile, probs=c(0.95, 0.05), R=600, na.rm=TRUE)
lines(quant.means[1,], col="blue", lty=2)
lines(quant.means[2,], col="blue", lty=2)
#plot mean and CI for historic songs.
#plot original song values
library("boot")
sample.mean <- function(d, x) {
mean(d[x])
}
boot_hist <- boot(starting.trait, statistic=sample.mean, R=100)#, strata=mn.res$iteration)
ci.hist <- boot.ci(boot_hist, conf=0.95, type="basic")
low <- ci.hist$basic[4]
high <- ci.hist$basic[5]
points(0, mean(starting.trait), pch=20, cex=0.6, col="black")
library("Hmisc")
errbar(x=0, y=mean(starting.trait), high, low, add=TRUE)
#plot current song values
boot_curr <- boot(target.data, statistic=sample.mean, R=100)#, strata=mn.res$iteration)
ci.curr <- boot.ci(boot_curr, conf=0.95, type="basic")
low <- ci.curr$basic[4]
high <- ci.curr$basic[5]
points(years, mean(target.data), pch=20, cex=0.6, col="black")
errbar(x=years, y=mean(target.data), high, low, add=TRUE)
#plot panel title
text(x=3, y=max(par.sens1$sens.results[ , , "trait.pop.mean", ], na.rm=TRUE)-100,
labels=paste("Par = ", par.range[i], sep=""))
}
SongEvo documentation built on May 1, 2019, 11:30 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value See Also Examples
Description
This function follows par.sens to help users optimize values for imperfectly known parameters for SongEvo. The goals are to maximize accuracy and precision of model prediction.
Usage
1 |
Arguments
sens.results |
The sens.results array from par.sens(), which includes summary.results from SongEvo() for a range of parameter values. summary.results from SongEvo() includes population summary values for each time step (dimension 1) in each iteration (dimension 2) of the model. Population summary values are contained in five additional dimensions: population size for each time step of each iteration (<e2><80><9c>sample.n<e2><80><9d>), the population mean and variance of the song feature studied (<e2><80><9c>trait.pop.mean<e2><80><9d> and <e2><80><9c>trait.pop.variance<e2><80><9d>), with associated lower (<e2><80><9c>lci<e2><80><9d>) and upper (<e2><80><9c>uci<e2><80><9d>) confidence intervals. |
ts |
The timestep (<e2><80><9c>ts<e2><80><9d>) at which to compare simulated trait values to target trait values (<e2><80><9c>target.data<e2><80><9d>). |
target.data |
A set of trait values against which to compare simulated values. target.data may be measured (e.g. from a training population) or hypothetical. |
par.range |
Range of parameter values over which to optimize. |
Value
Three measurements of accuracy and one measure of precision. Accuracy is quantified by three different approaches: i) the mean of absolute residuals of the predicted population mean values in relation to observed values (smaller absolute residuals indicate a more accurate model), ii) the difference between the bootstrapped mean of predicted population means and the mean of the observed values, and iii) the proportion of simulated population trait means that fall within confidence intervals of the observed data (a higher proportion indicates greater accuracy). Precision is measured with the residuals of the predicted population variance to the variance of observed values (smaller residuals indicate a more precise model).
See Also
SongEvo, par.sens, mod.val, h.test
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 | ### See vignette for an example that uses all functions in SongEvo.
#### Specify and call `par.sens()`
# Here we test the sensitivity of the Acquire a Territory submodel to variation
# in territory turnover rates, ranging from 0.8<e2><80><93>1.2 times the published rate
# (40<e2><80><93>60% of territories turned over). The call for the par.sens function has a
# format similar to SongEvo. The user specifies the parameter to test and the
# range of values for that parameter. The function currently allows examination
# of only one parameter at a time and requires at least two iterations.
parm <- "terr.turnover"
par.range = seq(from=0.45, to=0.55, by=0.05)
sens.results <- NULL
data("song.data")
data("glo.parms")
years=2005-1969
iteration=5
timestep=1
n.territories <- glo.parms$n.territories
starting.trait <- subset(song.data, Population=="PRBO" & Year==1969)$Trill.FBW
starting.trait2 <- c(starting.trait, rnorm(n.territories-length(starting.trait),
mean=mean(starting.trait), sd=sd(starting.trait)))
init.inds <- data.frame(id = seq(1:n.territories), age = 2, trait = starting.trait2)
init.inds$x1 <- round(runif(n.territories, min=-122.481858, max=-122.447270), digits=8)
init.inds$y1 <- round(runif(n.territories, min=37.787768, max=37.805645), digits=8)
# Now we call the par.sens function with our specifications.
extra_parms <- list(init.inds = init.inds,
timestep = 1,
n.territories = nrow(init.inds),
learning.method = "integrate",
integrate.dist = 0.1,
lifespan = NA,
terr.turnover = 0.5,
mate.comp = FALSE,
prin = FALSE,
all = TRUE)
global_parms_key <- which(!names(glo.parms) %in% names(extra_parms))
extra_parms[names(glo.parms[global_parms_key])]=glo.parms[global_parms_key]
par.sens1 <- par.sens(parm = parm, par.range = par.range,
iteration = iteration, steps = years, mate.comp = FALSE,
fixed_parms=extra_parms[names(extra_parms)!=parm], all = TRUE)
#### Prepare current song values
target.data <- subset(song.data, Population=="PRBO" & Year==2005)$Trill.FBW
#### Specify and call `par.opt()`
# Users specify the timestep (<e2><80><9c>ts<e2><80><9d>) at which to compare simulated trait values
# to target trait data (<e2><80><9c>target.data<e2><80><9d>) and save the results in an object (called
# `par.opt1` here).
ts <- years
par.opt1 <- par.opt(sens.results=par.sens1$sens.results, ts=ts,
target.data=target.data, par.range=par.range)
# Examine results objects (residuals and target match).
par.opt1$Residuals
par.opt1$Target.match
#### Plot results of `par.opt()`
#### Accuracy
#1. Difference in means.
plot(par.range, par.opt1$Target.match[,1], type="l",
xlab="Parameter range", ylab="Difference in means (Hz)")
#2. Plot proportion contained.
plot(par.range, par.opt1$Prop.contained, type="l",
xlab="Parameter range", ylab="Proportion contained")
#3. Calculate and plot mean and quantiles of residuals of mean trait values.
res.mean.means <- apply(par.opt1$Residuals[, , 1], MARGIN=1,
mean, na.rm=TRUE)
res.mean.quants <- apply (par.opt1$Residuals[, , 1], MARGIN=1,
quantile, probs=c(0.975, 0.025), R=600, na.rm=TRUE)
plot(par.range, res.mean.means, col="orange",
ylim=range(par.opt1$Residuals[,,1], na.rm=TRUE),
type="b",
xlab="Parameter value (territory turnover rate)",
ylab="Residual of trait mean (trill bandwidth, Hz)")
points(par.range, res.mean.quants[1,], col="orange")
points(par.range, res.mean.quants[2,], col="orange")
lines(par.range, res.mean.quants[1,], col="orange", lty=2)
lines(par.range, res.mean.quants[2,], col="orange", lty=2)
#### Precision
#Calculate and plot mean and quantiles of residuals of variance of trait values
res.var.mean <- apply(par.opt1$Residuals[, , 2], MARGIN=1,
mean, na.rm=TRUE)
res.var.quants <- apply (par.opt1$Residuals[, , 2], MARGIN=1,
quantile, probs=c(0.975, 0.025), R=600, na.rm=TRUE)
plot(par.range, res.var.mean, col="purple",
ylim=range(par.opt1$Residuals[,,2], na.rm=TRUE),
type="b",
xlab="Parameter value (territory turnover rate)",
ylab="Residual of trait variance (trill bandwidth, Hz)")
points(par.range, res.var.quants[1,], col="purple")
points(par.range, res.var.quants[2,], col="purple")
lines(par.range, res.var.quants[1,], col="purple", lty=2)
lines(par.range, res.var.quants[2,], col="purple", lty=2)
#### Visual inspection of accuracy and precision: plot trait values for range of parameters
par(mfcol=c(3,2),
mar=c(4.1, 4.1, 1, 1),
cex=1.2)
for(i in 1:length(par.range)){
plot(par.sens1$sens.results[ , , "trait.pop.mean", ],
xlab="Year", ylab="Bandwidth (Hz)",
xaxt="n", type="n",
xlim=c(-0.5, years), ylim=range(par.sens1$sens.results[ , , "trait.pop.mean", ], na.rm=TRUE))
for(p in 1:iteration){
lines(par.sens1$sens.results[p, , "trait.pop.mean", i], col="light gray")
}
freq.mean <- apply(par.sens1$sens.results[, , "trait.pop.mean", i], 2, mean, na.rm=TRUE)
lines(freq.mean, col="blue")
axis(side=1, at=seq(0, 35, by=5), labels=seq(1970, 2005, by=5))#, tcl=-0.25, mgp=c(2,0.5,0))
#Plot 95% quantiles
quant.means <- apply (par.sens1$sens.results[, , "trait.pop.mean", i], MARGIN=2,
quantile, probs=c(0.95, 0.05), R=600, na.rm=TRUE)
lines(quant.means[1,], col="blue", lty=2)
lines(quant.means[2,], col="blue", lty=2)
#plot mean and CI for historic songs.
#plot original song values
library("boot")
sample.mean <- function(d, x) {
mean(d[x])
}
boot_hist <- boot(starting.trait, statistic=sample.mean, R=100)#, strata=mn.res$iteration)
ci.hist <- boot.ci(boot_hist, conf=0.95, type="basic")
low <- ci.hist$basic[4]
high <- ci.hist$basic[5]
points(0, mean(starting.trait), pch=20, cex=0.6, col="black")
library("Hmisc")
errbar(x=0, y=mean(starting.trait), high, low, add=TRUE)
#plot current song values
boot_curr <- boot(target.data, statistic=sample.mean, R=100)#, strata=mn.res$iteration)
ci.curr <- boot.ci(boot_curr, conf=0.95, type="basic")
low <- ci.curr$basic[4]
high <- ci.curr$basic[5]
points(years, mean(target.data), pch=20, cex=0.6, col="black")
errbar(x=years, y=mean(target.data), high, low, add=TRUE)
#plot panel title
text(x=3, y=max(par.sens1$sens.results[ , , "trait.pop.mean", ], na.rm=TRUE)-100,
labels=paste("Par = ", par.range[i], sep=""))
}
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.