Description Usage Arguments Examples
This function uses samples from a stanfit object to generate predictive checks on directions of travel.
1 2 3 | postPredAngle(model.fit, df.obs = NULL, contours = c(25, 50, 75),
rangePred = 1:10, numbDraws = 500, buffer = 10, extention = 10,
angles = NULL)
|
model.fit |
stanfit model, containing a_pred as generated quantities: a_pred is the model predicted angle. |
df.obs |
dataframe of observed travel points. |
contours |
Contours used in the plot of the kernel density of predicted points, if null a raster is produced. |
rangePred |
Range of observed travel points to predict. |
numbDraws |
Number of predictions to make for each observed point. |
buffer |
Display parameter, used to extend the plot extent range. |
extention |
Display parameter, used to define the radius on which to plot predictions for each point. |
context |
List of shapefiles or rasters that can provide contextual information. |
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 | ## Not run:
library(rstan)
#################################
# Get data ready for rStan #
#################################
focalBaboon[is.na(focalBaboon)]<-0
data.for.stan<-list(N = nrow(focalBaboon),
ax = focalBaboon$dx.obs,
ay = focalBaboon$dy.obs,
bx=focalBaboon$dx.bearing,
by=focalBaboon$dy.bearing,
ax_cv=focalBaboon$dx.resultant,
ay_cv=focalBaboon$dy.resultant,
ax0=focalBaboon$dx.0,
ax1=focalBaboon$dx.1,
ay0=focalBaboon$dy.0,
ay1=focalBaboon$dy.1)
#################################
# Define rStan model #
#################################
stan.model.test = '
data{
int<lower = 0> N; //number of data points
vector[N] ax; //displacement of observed travel
vector[N] ay; //displacement of observed travel
vector[N] bx; //displacement of previous travel
vector[N] by; //displacement of previous travel
vector[N] ax_cv; //X component of unit vector towards group
vector[N] ay_cv; //Y component of unit vector towards group
vector[N] ax0; //X component of unit vector towards individual 0
vector[N] ay0; //Y component of unit vector towards individual 0
vector[N] ax1; //X component of unit vector towards individual 1
vector[N] ay1; //Y component of unit vector towards individual 1
}
transformed data{
vector[N] a;
//set observed angle of travel
for(i in 1:N){
a[i] = atan2(ay[i],ax[i]);
}
}
parameters{
//These are the parameters to be estimated for angle of travel
real<lower = 0> kappa;
real beta_cv;
real beta_a0;
real beta_a1;
}
model{
vector[N] y_hat;
//setting our priors for betas
kappa ~ normal(0,10);
beta_cv ~ normal(0,1);
beta_a0 ~ normal(0,1);
beta_a1 ~ normal(0,1);
//running the model
for(i in 1:N){
//estimating angle of travel
y_hat[i] = atan2( by[i] + beta_cv*ay_cv[i] + beta_a0*ay0[i] + beta_a1*ay1[i], bx[i] + beta_cv*ax_cv[i] + beta_a0*ax0[i] + beta_a1*ax1[i]);
}
//likelihood
a ~ von_mises(y_hat, kappa);
}
generated quantities{
//generate log-likelihood for observations & predictions
vector[N] log_lik;
vector[N] a_pred;
for(i in 1:N){
a_pred[i] = von_mises_rng(atan2( by[i] + beta_cv*ay_cv[i], bx[i] + beta_cv*ax_cv[i] ),kappa);
log_lik[i] = von_mises_lpdf(a[i]|atan2( by[i] + beta_cv*ay_cv[i], bx[i] + beta_cv*ax_cv[i] ),kappa);
}
}
'
#################################
# Fit the model #
#################################
fit.1 = stan(model_code = stan.model.test, data = data.for.stan, iter=1000, chains=1, cores=1)
#################################
# Check model fit #
#################################
print(fit.1, digits = 4, pars=c("kappa","beta_cv","beta_a0","beta_a1"))
#launch_shinystan(fit.1)
#################################
# Posterior predictive checks #
#################################
postPredAngle(fit.1, focalBaboon, rangePred=1:10)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.