Last updated: 2018-05-04
Code version: fc4400b
We will use R version 3.3.2 for compatibility
library(dplyr)
library(data.table)
library(magrittr)
library(purrr)
library(here) # for tracking working directory
library(ggplot2)
library(epistats)
library(broom)
library(rjags)
library(Bain)
Several studies suggest that cognitive behavioral therapy is an effective treatment for Posttrau- matic Stress Disorder (PTSD) in male veterans. Imagine that you performed a study to compare Prolonged Exposure (PE), a type of cognitive behavioral therapy, with Present-Centered therapy (PC), a supportive intervention. In this randomized controlled trial 284 veterans were assigned to receive either PE or PC. The (dichotomous) outcome measure of interest was “loss of diagnosis”(LD); that is, each veteran that was no longer diagnosed with PTSD after treatment got the score 1, and the ones that still had the diagnosis got the score 0.
In this exercise you will use rjags to examine whether PE and PC are effective treatments for PTSD in veterans. The data you need for this exercise are displayed in Table 1 below.
Table 1 PE PC Loss of diagnosis 58 40 Total men treated 141 143
Open the file Model1.txt (you can open it in R, or in a text editor like notepad). We will specify the model in this file in a couple of steps. In the the Model1.txt file you can add your own comments by putting # in front of those comments. In that way, JAGS will know that it should not run that code, but consider it as a comment
We are interested in the proportion of veterans that had a LD in each treatment condition, that is, we are interested in the probability of observing a succes (a 1) in the pe condition (θ_pe) and the pc condition (θ_pc).
We can model the number of successes for each condition with binomial distributions:
y_pe ~ Binomial(θ_pe, n_pe)
y_pc ~ Binomial(θ_pc, n_pc)
That is, for each condition, we assume that the data y (the number of successes in the data) have a binomial distribution with as parameters θ, the probability of observing a success, and n, the total number of observations in the condition.
The total number of observations in the conditions we know, as well as of course the data y (see Table 1). The probabilities θ¸ are what we want to estimate from the information we have, and we will do this in a Bayesian way.
In the model code for JAGS in Model1.txt we have to specify these distributions within model{ } so that JAGS can find the likelihood based on this information. Therefore, within model{ } we add the following code in Model1.txt:
likelihood of the data
ype ~ dbin(ppe,npe) ypc ~ dbin(ppc,npc) Step 2. Priors Next, we have to specify prior distributions for all the parameters we need to estimate. That is, we need priors for θ_pe (ppe in our code) and θ_pc (ppc in our code). A commonly used conjugate prior distribution for probabilities from a binomial distribution is the beta-distribution. That is,
θ_pe ~ Beta(a_pe, b_pe)
θ_pc ~ Beta(a_pc, b_pc)
with shape parameters a (a_pe, a_pc) and b (b_pe, b_pc).
For each condition, we will specify that a = b = 1 in our model code to obtain low-informative prior distributions. That is, we add the following code in Model1.txt within model{ }:
prior distributions
ppe ~ dbeta(1,1) ppc ~ dbeta(1,1) Step 3. Interesting additional quantities Finally we will specify a contrast you may be interested in, namely the ratio of the probability of succes in condition PE and the probability of succes in condition PC (θ_pe /θ_pc). By adding the following code within model{ } we will tell JAGS to calculate this ratio in each iteration of the estimator:
contrast
ratio <- ppe/ppc You could also add any other quantity based on the probabilities in your code, for example, the difference in the probabilities.
Our model code is done for now - save your code as Model1.txt in your working directory.
Before we can use our model, we first have to get our data (from Table 1) into R. Make four variables in R that contain respectively the number of veterans that lost their diagnosis for condition PE, the number of veterans that lost their diagnosis for condition PC, the total number of veterans that were in condition PE, and the total number of veterans that were in condition PC.
data for exercise 1
y_pe <- 58
y_pc <- 40
n_pe <- 141
n_pc <- 143
Now, if we want to use this data for our analysis in rjags we need to put these data in what is called a list. In the list we specify for all the names we used for the data in the jags model code (that is, ype, ypc, npe, npc) what the actual data is that we have stored in R (y_pe, y_pc, n_pe, n_pc).
About lists: You can store various types of R objects (such as vectors, matrices or other lists) in a list. The list type of object is used for specifying data and initial values in rjags amongst other things, but it is, for instance, also frequently used for storing and presenting output from analyses in R. If you want to know more about lists you can read more here: http://rforpublichealth.blogspot.nl/2015/03/basics-of-lists.html.
Make a list out of the data, and call the list you made to see what is inside as follows:
data for exercise 1 for rjags
(actually, naming arguments in lists works in R without parentheses)
dat_ex1 <- list('ype'=y_pe, 'ypc'=y_pc, 'npe'=n_pe, 'npc'=n_pc)
dat_ex1
$ype
[1] 58
$ypc
[1] 40
$npe
[1] 141
$npc
[1] 143
dat_ex1 <- list(ype=y_pe, ypc=y_pc, npe=n_pe, npc=n_pc)
dat_ex1
$ype
[1] 58
$ypc
[1] 40
$npe
[1] 141
$npc
[1] 143
In order to run the analyses we first have to tell rjags what the model file is, what the data is, and specify the number of chains for the analyses (we use 2 for this exercise - you will learn more about this on Tuesday). It is also possible to specify initial values, but we will skip this for now (JAGS will generate the initial values itself based on the priors you specified in Model1.txt). We specify the model by means of the function jags.model(). To learn a bit more about this function already, run ?jags.model in R.
We use the jags.model() function and we name it rjags_ex1 as follows :
nsamp = 10000
rjags_ex1 <- jags.model(file=here("analysis", "bayes_1_Model1.txt"),
data = dat_ex1,
n.chains = 2) ###compile and initialize jags model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2
Unobserved stochastic nodes: 2
Total graph size: 9
Initializing model
After you have run this line of code, jags will compile the model and specify initial values for (initialize) the model. Now, we want to start the estimation procedure to fit our Bayesian model. We start by letting the estimation procedure run for a while, to start up and converge (these first few iterations in the estimation procedure are called burn-in - more on this on Tuesday). We do this with the function update(). Using update(), jags starts the estimation procedure, but it doesn’t save this as the results.
We tell update which model to update in the first argument, and how many times in the second argument as follows:
update(rjags_ex1, nsamp) ####burnin
When jags is done with the first updates, we can now tell it to run some more iterations in our estimation procedure, and to calculate and store some results for us. For this we use the function coda.samples(). We need to tell coda.samples() which model to update, for which parameters we want jags to save the results, and how many iterations to use (we use 1000). To learn more about this function, use ?coda.samples().
First we will make a vector with function c() with the names of all the parameters we want to see the results for. Next, we specify the function coda.samples, and we name it samples_rjags_ex1.
parstosave=c('ppe','ppc', 'ratio')
samples_rjags_ex1=coda.samples(model=rjags_ex1, variable.names=parstosave, n.iter=nsamp)
When it has reached 100% it is done. For this model, it will go really fast, but for more complicated models it might take a little while, depending on how many iterations you chose.
After updating you normally have to check whether the algorithm has reached the target distribution. On Tuesday you will learn why and how this is done. For now, you may assume the sampler has reached convergence, and you can go on with analyzing the posterior results.
If you call samples_rjags_ex1 you will get all the iterations for all the parameters rjags has saved for you. It can be very convenient to have all of these samples, however, for now we just want to see some summary results. To get summary results use function summary() on samples_rjags_ex1.
summary(samples_rjags_ex1)
Iterations = 11001:21000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
ppc 0.2829 0.03720 0.0002631 0.0003266
ppe 0.4123 0.04099 0.0002898 0.0003654
ratio 1.4835 0.24980 0.0017663 0.0022190
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
ppc 0.2129 0.2571 0.2820 0.3081 0.3579
ppe 0.3338 0.3843 0.4115 0.4403 0.4928
ratio 1.0634 1.3078 1.4589 1.6356 2.0371
You will see the mean, standard deviation, and various quantiles of the posteriors for each estimated parameter (you will also see two kinds of standard errors, but you can disregard these for now). These summary statistics are calculated by rjags over all of the saved samples per parameter in samples_rjags_ex1.
You can also makes plots of the posterior distributions like this:
plot(samples_rjags_ex1,density=TRUE, trace=FALSE) ###densities
Write down and interpret these posterior results (mean, median, 95% credibility interval). Do you think PE is an effective treatment for PTSD in male veterans? What about PC? Would you consider one treatment superior to the other? What do you base this on?
Based on the ratio of the two, the 95% credible interval does not include the null-value of no effect (ratio = 1), so it seems that ppe is more effective
Rjags has provided you with the posterior results for θ_pe, θ_pc, and their ratio. In this part of the exercise you will obtain the posterior means for these parameters analytically.
It can be shown that the posterior distribution for each θ i is given by: It can further be shown that the posterior mean ( θ̂ i ) of the distribution is given by:
Calculate the posterior mean for both proportions θ_pe and θ_pc using this equation. Then calculate the ratio of these means. Compare these results to those you have obtained using rjags. Are the results similar?
a = 1; b = 1
theta_pe <- (a + y_pe) / (a + y_pe + b + n_pe - y_pe)
theta_pc <- (a + y_pc) / (a + y_pc + b + n_pc - y_pc)
theta_pe
[1] 0.4125874
theta_pc
[1] 0.2827586
theta_pe / theta_pc
[1] 1.459151
On an international conference on PTSD you meet two fellow researchers, Thomas B. and Ronald F., who also evaluated the use of PE versus PC. Thomas recently executed a randomized clinical trial with 520 female veterans to evaluate the use of PE versus PC. Ronald F. did a comparable trial with 235 WOII male veterans in 1946.
As discussed during the today’s lecture on informative prior specification, it can we worthwhile to include data obtained in previous studies in the analysis of new data.
Would you be interested to include the results obtained in either of the trials in the prior distribution for the analysis of your own data? Which data do you think would be most relevant and why?
thomas has data on only women, which is a small subset of all veterans, and women differ from men with regards to psychiologic disorders
ronald used data from 1946, which is a very long time ago
I would give more weight to the ronald’s data
so use these beta distributions: ppe ~ dbeta(41,66) ppc ~ dbeta(45,86)
So now:
40 / 105
[1] 0.3809524
45 / 130
[1] 0.3461538
(40 / 105) / (45 / 130)
[1] 1.100529
In Exercise 2.c you will rerun the analysis from Exercise 1 with informative prior distributions based on the data obtained by either Thomas or Ronald. In order to do so you will need to alter the uninformative Beta(1,1) distributions into informative distributions. This is what their data looks like:
Table 2. PTSD data from Thomas B. PE PC Loss of diagnosis 120 80 Total men treated 245 275 Table 3. PTSD data from Ronald F. PE PC Loss of diagnosis 40 45 Total men treated 105 130 Given the data that you consider to be the most relevant, what would your new informative prior distributions look like?
Hint: You can think of parameters a_i and b_i of the Beta-distribution as the prior number of successes and failures + 1 respectively.
Run the analysis again with your informative prior distributions. Ask for the posterior results for both proportions θ_pe and θ_pc and their ratio. Compare the results (mean, median, 95% Central Credibility Interval) with the results obtained in Exercise 1. To what degree do the informative priors influence the posterior results? Do you think this is an desirable effect?
rjags_ex2 <- jags.model(file=here("analysis", "bayes_1_Model2.txt"),
data = dat_ex1,
n.chains = 2) ###compile and initialize jags model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2
Unobserved stochastic nodes: 2
Total graph size: 12
Initializing model
update(rjags_ex2, nsamp)
samples_rjags_ex2 = coda.samples(model=rjags_ex2, variable.names=parstosave, n.iter=nsamp)
summary(samples_rjags_ex2)
Iterations = 11001:21000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
ppc 0.3128 0.0280 0.0001980 0.0002526
ppe 0.3992 0.0309 0.0002185 0.0002913
ratio 1.2867 0.1531 0.0010828 0.0013985
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
ppc 0.2596 0.2938 0.3121 0.3312 0.3698
ppe 0.3401 0.3782 0.3989 0.4201 0.4601
ratio 1.0113 1.1798 1.2784 1.3828 1.6128
The ratio is lower, but the 95% credible interval still does not include 1
Rerun the analysis with informative prior distributions based on the data obtained with the other randomized controlled trial. Compare the results from all three analyses. Which informative priors affect the results the most? Why?
120 / 245
[1] 0.4897959
80 / 275
[1] 0.2909091
(120 / 245) / (80 / 275)
[1] 1.683673
rjags_ex2.2 <- jags.model(file=here("analysis", "bayes_1_Model2.2.txt"),
data = dat_ex1,
n.chains = 2) ###compile and initialize jags model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2
Unobserved stochastic nodes: 2
Total graph size: 12
Initializing model
update(rjags_ex2.2, nsamp)
samples_rjags_ex2.2 = coda.samples(model=rjags_ex2.2, variable.names=parstosave, n.iter=nsamp)
summary(samples_rjags_ex2.2)
Iterations = 11001:21000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
ppc 0.2863 0.02220 0.0001570 0.0001939
ppe 0.4613 0.02568 0.0001816 0.0002337
ratio 1.6212 0.15658 0.0011072 0.0013423
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
ppc 0.2441 0.2709 0.2860 0.3011 0.3310
ppe 0.4109 0.4438 0.4613 0.4784 0.5118
ratio 1.3379 1.5132 1.6124 1.7217 1.9515
Using these priors, the ratio comes out higher
summaries <- map(list(model_1 = samples_rjags_ex1,
model_2 = samples_rjags_ex2,
model_2.2 = samples_rjags_ex2.2),
summary)
map(summaries, "quantiles")
$model_1
2.5% 25% 50% 75% 97.5%
ppc 0.2128730 0.2570714 0.2819624 0.3080776 0.3578587
ppe 0.3338448 0.3843227 0.4114886 0.4403394 0.4927564
ratio 1.0633575 1.3078245 1.4588771 1.6355836 2.0370514
$model_2
2.5% 25% 50% 75% 97.5%
ppc 0.2596059 0.2937855 0.3120627 0.3312302 0.3698253
ppe 0.3400957 0.3781606 0.3988612 0.4201098 0.4600784
ratio 1.0112539 1.1798337 1.2783770 1.3827924 1.6127518
$model_2.2
2.5% 25% 50% 75% 97.5%
ppc 0.2441026 0.2709388 0.2859529 0.3010749 0.3309539
ppe 0.4109234 0.4437891 0.4613238 0.4784141 0.5117843
ratio 1.3378968 1.5132265 1.6124217 1.7216934 1.9514837
The data from Thomas were based on more people, but whe estimated coefficient was closer to what we observed in our data, so the influence on the estimate was smaller
What happens when you include conflicting data? What happens when your prior is based on more data than available in the current study?
It overrides the data
Remember we used ppe ~ dbeta(41,66) ppc ~ dbeta(45,86)
And
ppe ~ dbeta(121,126) ppc ~ dbeta(80,196)
a_pe = 41; b_pe = 66
a_pc = 46; b_pc = 86
theta_pe <- (a_pe + y_pe) / (a_pe + y_pe + b_pe + n_pe - y_pe)
theta_pc <- (a_pc + y_pc) / (a_pc + y_pc + b_pc + n_pc - y_pc)
theta_pe / theta_pc
[1] 1.276491
theta_pe
[1] 0.3991935
theta_pc
[1] 0.3127273
a_pe = 121; b_pe = 126
a_pc = 80; b_pc = 196
theta_pe <- (a_pe + y_pe) / (a_pe + y_pe + b_pe + n_pe - y_pe)
theta_pc <- (a_pc + y_pc) / (a_pc + y_pc + b_pc + n_pc - y_pc)
theta_pe / theta_pc
[1] 1.610846
How did you include the prior information in the calculation of the posterior mean? Did you get similar results?
Yes the results are very similar.
In fact, they were closer to the 50% percentile than to the means
map(summaries, "statistics")
$model_1
Mean SD Naive SE Time-series SE
ppc 0.2828738 0.03720338 0.0002630677 0.0003266170
ppe 0.4123127 0.04098871 0.0002898339 0.0003653588
ratio 1.4835232 0.24979780 0.0017663372 0.0022190393
$model_2
Mean SD Naive SE Time-series SE
ppc 0.3127617 0.02799939 0.0001979856 0.0002525555
ppe 0.3992100 0.03089787 0.0002184810 0.0002912842
ratio 1.2866978 0.15313240 0.0010828096 0.0013984623
$model_2.2
Mean SD Naive SE Time-series SE
ppc 0.2862550 0.02220056 0.0001569817 0.0001938545
ppe 0.4612668 0.02568016 0.0001815862 0.0002336660
ratio 1.6212210 0.15657696 0.0011071663 0.0013423318
In this exercise you will focus on the data from the PE-group. In the previous exercise you assumed these data were binomially distributed. In this exercise you will assess whether the binomial model actually fits the data. If all the observed statuses of diagnosis are truly a sequence of independent Bernoulli trials (such that the number of successes is binomially distributed), the proportion of success in the first half of the data should be equal to the proportion of success in the second half. You will evaluate this using rjags, by performing a posterior predictive check.
Specify the basic model for the bernouilli distributed data of the PE condition, with uninformative prior distributions.
Open the file Model2.txt. Within the model{} statement, we need to specify the distributions and equations that describe the data for the PE-condition, so that JAGS can determine the likelihood of the data. The data you will use for this exercise is a series of n=141 Bernoulli trials. Please refer to the data in your model with x (like how in the previous exercises the data was called y_pe or y_pc), and refer to the total sample size with n.
You can do this with the following code: for (i in 1:n){ x[i] ~ dbern(ppe) } You see this code makes use of a “for loop”. Inside the for loop, it says x[i] ~ dbern(ppe). The for loop works like this: In its first loop, it note that i is equal to 1. Then it goes to the equation x[i] ~ dbern(ppe), and reads this as x[1] ~ dbern(ppe), that is, the first observation in the data x comes from a bernoulli distribution with parameter ppe. Then, it finishes its first loop, and goes to the second loop. Now it notes that i=2, and that x[2] ~ dbern(ppe), that is, the second observation in the data x also comes from a bernoulli distribution with parameter ppe. Then it goes to the third loop, and so on, until i is equal to n (with n being equal to the sample size of x). That is, this code says that each observation in x comes from a bernoulli distribution with parameter ppe.
Now, in file Model2.txt, within the model{} statement, specify a Beta(1,1) prior distribution for θ_pe.
Specify all additional quantities in your model file that you need to do the posterior predictive check.
Obtain the observed proportion of successes in each half of your data, and the difference between these proportions Within the model{} statement specify you wish to calculate the proportion of success in the first and second half of the data. Then specify you wish to calculate the difference between the two proportions. The latter is your discrepancy measure of interest.
You can do this with the following code:
discrepancy measure in the data
ppe1 <- sum(x[1:70])/70 ppe2 <- sum(x[71:141])/71 dif <- ppe1-ppe2 More information on the sum-function and other logical functions can be found in the JAGS manual under “Functions”.
Obtain simulated data sets that are in line with your model assumptions. By sampling from a binomial distribution with the same probability as in the real dataset, we can obtain simulated data sets of bernoulli trials (each sample from the binomial distribution is equal to one data set of bernoulli trials). Specify two binomial distributions - the first one with a sample size that is equal to the first half of that of the real data, the second one with the same size as the second half of the real data. If the real data is also truly a set of bernouilli trials, the proportions of succes in the simulated datasets should be similar to what we found for the real dataset.
You can thus add the following lines within the model{} statement to obtain samples from the posterior predictive distribution in each iteration of the sampler (replicates of the first and second halves of the data):
posterior predictive distribution
postpred.ype1 ~ dbin(ppe,70) postpred.ype2 ~ dbin(ppe,71) Step 3. Obtain the observed proportion of successes in each half of each simulated data set, and the difference between these proportions. Obtain the proportion of success in the first and second half of the replicated data (simular to step 1), and calculate the difference between the two proportions (you could call it postpred.dif). This piece of code is not given, so you have to program this yourself.
Step 4. Use the step function to obtain the posterior predictive p-value (ppp-value) You can use the step function in JAGS in the following way to get a posterior predictive p-value for the difference between the real and simulated difference between the proportions of successes in each half of the data:
posterior predictive p-value
p <- step(postpred.dif - dif) What does the step function do? See the JAGS manual (in chapter ‘Functions’) or (Ntzoufras 2009, p.96) for details on the step() function.
Step function is False when x < 0, and True when x >= 1
Now you know what the step function does, what does the mean of the posterior of p mean based on this piece of code? What is H0 in our (posterior predictive) hypothesis? What would you expect the posterior distribution of p to look like if H0 were true?
In what proportion of the simulated samples was the difference higher than the observed difference.
If H0 is true, the observed difference should fall within the bulk of the distribution of the simulated differences
Run this line of code in R to obtain a list with the data for the rjags model:
dat_ex3= list(x=c(0,0,1,0,1,0,1,1,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,0,0,1,0,1,1,0,0,0,1,0,1,0,1,1,1,1,1,0,0,1,1,1,0,1,0,1,1,0,0,0,1,0,0,1,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,1,0,0,0,0,0,1,0,0,1,1,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,0,1,1,0,0,0,1,0,0,0,0,0,0,1,0,1,0,0,0,1,1,0,0,0,1,0,0,1,0,1,0,0,0,0,1,1,0,0)
, n=141)
Now run the rjags analysis on your data with the model you specified in Model2.txt, in the same general way as in Exercise 1 and 2.
nsamp = 1000
rjag_3 <- jags.model(here("analysis", "bayes_1_Model3.txt"),
data = dat_ex3,
n.chains = 2)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 141
Unobserved stochastic nodes: 3
Total graph size: 162
Initializing model
update(rjag_3, nsamp)
vars_to_save <- c("postpred.ppe1", "postpred.ppe2", "postpred.dif", "p")
rjag_3_samples <- coda.samples(rjag_3, variable.names = vars_to_save, n.iter = nsamp)
summary(rjag_3_samples)
Iterations = 1001:2000
Thinning interval = 1
Number of chains = 2
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
p 0.214500 0.41058 0.009181 0.009180
postpred.dif -0.003471 0.08129 0.001818 0.001788
postpred.ppe1 0.410564 0.07181 0.001606 0.001606
postpred.ppe2 0.414035 0.07057 0.001578 0.001578
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
p 0.0000 0.00000 0.000000 0.00000 1.0000
postpred.dif -0.1642 -0.06358 -0.007847 0.04889 0.1620
postpred.ppe1 0.2714 0.35714 0.414286 0.45714 0.5571
postpred.ppe2 0.2817 0.36620 0.408451 0.46479 0.5634
What are your conclusions? Finally, give some thought to how this relates to the classical alternative for the ppp-value. Describe the difference(s).
The proportion of times the posterior predictive check difference was higher than the observed difference
This was in 22% of the cases, which means that the observed difference is not very extreme under the null-hypothesis
x1 <- sum(dat_ex3$x[1:70])
n1 <- 70
x2 <- sum(dat_ex3$x[71:141])
n2 <- 71
prop.test(c(x1, x2), c(n1, n2))
2-sample test for equality of proportions with continuity
correction
data: c(x1, x2) out of c(n1, n2)
X-squared = 0.34087, df = 1, p-value = 0.5593
alternative hypothesis: two.sided
95 percent confidence interval:
-0.1137550 0.2389059
sample estimates:
prop 1 prop 2
0.4428571 0.3802817
The p-values are not very close
Considerations for your presentation on Friday: Perform a Bayesian analysis on some (your own?) data set with uninformative and informative prior distributions Perform a posterior predictive pheck to check model assumptions for a particular (your own?) data set Present about the difference between classical hypothesis testing and using posterior predictive p-values
Applied Bayesian Statistics: Lab Meeting 2. Convergence, model selection with the DIC, and prediction. Today’s objectives Today you will use rjags to sample from the posterior distribution of a logistic regression model, and see if the procedure seems to have converged. You will also practice using the information criterium DIC in the context of the selection of predictors for a logistic regression model. Furthermore, for a person with known predictors but unknown outcome, you will obtain a prediction and credibility interval for the unknown outcome.
In this exercise you will develop a clinical prediction rule for detecting major depressive disorder in primary care. The dataset you will use for this exercise consists of various measures of 1046 participants, aged 18-65 years. These participants were recruited from general practice waiting rooms. Major depressive disorder (depr: 0=no depression, 1=depression) was assessed with the Composite International Diagnostic Interview according to DSM-IV-TR criteria. For this exercise we will predict depression diagnosis with gender (0=male, 1=female).
Open file ModelDay2.txt, and specify a univariate logistic regression model with gender as a predictor variable, such that JAGS can determine the likelihood of the data.
A logistic regression model could be specified for OpenBUGS as follows:
likelihood of the data
for (i in 1:n){ logit(p[i]) <- alpha + b.gender*gender[i] depr[i] ~ dbern(p[i])}
con1 <- file(here("analysis", "bayes_2_Modelday2.txt"), open = "r")
lines1 <- readLines(con1)
Next, specify prior distributions for all parameters in the model that need to be estimated. What type of prior distributions would you specify for the regression coefficient and for the intercept? How would you specify these distributions to be uninformative?
Take
You can refer to the jags manual for an overview of all available probability distributions (https://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/jags_user_manual.pdf/download). Pay special attention to what the parameters of each distribution represent - for example, the normal distribution is not specified with a mean and standard deviation as you might expect, but with a mean and a precision (the precision is equal to 1/variance).
You might want to specify the exponent of your regression coefficient in order to interpret directly the posterior distribution of the odds and odds ratios. To do so, include the following commands in the model:
odds ratios
odds0 <- exp(alpha) or.gender <- exp(b.gender)
With these lines of code you can load and specify the data for rjags:
load(here("data", "data_practical2.Rdata"))
data_prac2_ex1 <- dat_prac2[,1:2]
n=length(data_prac2_ex1[,1])
data_ex1 <- list("n"=n, "depr"= data_prac2_ex1[,1],"gender"= data_prac2_ex1[,2])
This time you will have also specify the initial values yourself (instead of generating initial values randomly). You should specify initial values for all the parameters in your model, and put them in a list. Make three of these lists, each with different initial values - one for each chain. Here is an example for the first of your chains:
Finally, put all three of your lists into a list like so:
ins1 <- list(alpha=0, b.gender=1)
ins2 <- list(alpha = -5, b.gender = 0)
ins3 <- list(alpha = 5, b.gender = -5)
ivals <- list(ins1, ins2, ins3)
Now specify your jags.model() function similar to the exercises from yesterday. Use three chains, and add inits= ivals as an argument to specify initial values for these chains.
jm_logit <- jags.model(here("analysis", "bayes_2_Modelday2.txt"),
data = data_ex1,
n.chains = length(ivals),
inits = ivals)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 2
Total graph size: 2105
Initializing model
Specify that you want to take 2000 samples burnin. Then specify which parameters you want to save, and use coda.samples() to take 2000 samples (per chain) to base your results on. Store the latter samples in object samps_ex1.
nburnin = 2000
update(jm_logit, nburnin)
nsamps = 2000
to_save <- c("alpha", "b.gender", "odds0", "or.gender")
samps_ex1 <- coda.samples(jm_logit, variable.names = to_save, n.iter = nsamps)
Inspect the trace plots You can obtain trace plots for each parameter with the following line of code:
plot(samps_ex1, density = F)
Do the three chains seem to mix properly?
Yes pretty good mix
Inspect the density plots You can obtain trace plots for each parameter with the following line of code:
plot(samps_ex1, trace = F)
Do the densities look nice and smooth? If not, you might want to take more samples.
Yes they are pretty smooth
Inspect the Gelman Rubin convergence diagnostics You can get both plots for the gelman rubin statistics for all the samples, and obtain summary statistics for them in the following way:
gelman.plot(samps_ex1) ##plot gelman rubin statistics
gelman.diag(samps_ex1) ## table gelman rubin statistics
Potential scale reduction factors:
Point est. Upper C.I.
alpha 1 1.01
b.gender 1 1.01
odds0 1 1.01
or.gender 1 1.01
Multivariate psrf
1
Does the black line in the plot converge to 1? Are the point estimates and majority of plausible values very close to 1 on equal to 1 for each parameter?
Yes all very close
Inspect the monte carlo standard errors (naive se) for each parameter Use the function summary() on samples_ex1 and inspect the standard errors for each parameter. The monte carlo standard errors (both the naive, and the time series error which is corrected for the autocorrelation of the chain) represent an estimate of the uncertainty contributed by only having a finite number of samples from the posterior. These errors should be small compared to the standard deviation of your estimates - a rule of thumb is that they should be smaller than about 5% of the posterior standard deviation. Is this the case for each parameter?
summary(samps_ex1)
Iterations = 3001:5000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha -2.0687 0.16316 0.0021063 0.0065371
b.gender 0.4867 0.19201 0.0024788 0.0078679
odds0 0.1280 0.02072 0.0002675 0.0008035
or.gender 1.6573 0.32446 0.0041888 0.0135047
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha -2.40952 -2.1729 -2.0633 -1.9568 -1.7576
b.gender 0.11785 0.3576 0.4815 0.6124 0.8791
odds0 0.08986 0.1138 0.1270 0.1413 0.1725
or.gender 1.12507 1.4299 1.6184 1.8448 2.4086
Yes they are lower than 5% from the SD of the parameter
Do you think the sampler has converged? If not, consider centering the predictor to avoid a strong correlation between the intercept and slope (as discussed in today’s lecture), and/or taking more samples.
To me it seems pretty converged
Tip: You can make a scatter plot to visually inspect the correlation between the intercept and slope for the first chain like this:
Make them into a data.frame, bind together by chain
samps_df <- samps_ex1 %>%
map(as.data.frame) %>%
rbindlist(idcol = "chain")
Plot by chain
samps_df %>%
ggplot(aes(x = alpha, y = b.gender)) +
geom_point(alpha = .5) +
facet_wrap(~chain)
Does not seem too narrow in either of the 3 chains
Once you have ensured yourself that the sampler has converged (for all parameters) you can continue to report and interpret the posterior results. What is the posterior estimate of the OR for gender? Looking at the 95% Credible Interval would you want to keep this predictor in the model?
Use the function summary() and your density plots to draw your conclusions.
summary(samps_ex1)
Iterations = 3001:5000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha -2.0687 0.16316 0.0021063 0.0065371
b.gender 0.4867 0.19201 0.0024788 0.0078679
odds0 0.1280 0.02072 0.0002675 0.0008035
or.gender 1.6573 0.32446 0.0041888 0.0135047
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha -2.40952 -2.1729 -2.0633 -1.9568 -1.7576
b.gender 0.11785 0.3576 0.4815 0.6124 0.8791
odds0 0.08986 0.1138 0.1270 0.1413 0.1725
or.gender 1.12507 1.4299 1.6184 1.8448 2.4086
From the 95% credible interval, I would keep the regression coefficient in.
Describe the difference between the classical Confidence Interval and the Bayesian Credibility Interval.
Check with glm
fit_glm <- glm(depr ~ gender, family = binomial, data = data_prac2_ex1)
summary(fit_glm)
Call:
glm(formula = depr ~ gender, family = binomial, data = data_prac2_ex1)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.6122 -0.6122 -0.6122 -0.4888 2.0899
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.0644 0.1638 -12.604 <2e-16 ***
gender 0.4850 0.1932 2.511 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 884.65 on 1045 degrees of freedom
Residual deviance: 878.02 on 1044 degrees of freedom
AIC: 882.02
Number of Fisher Scoring iterations: 4
confint(fit_glm)
2.5 % 97.5 %
(Intercept) -2.3996663 -1.7559256
gender 0.1138989 0.8728945
The confidence intervals are not the same.
Based on what you know about depressive disorders, select three candidate predictor variables (in addition to gender) from the variables listed below. Choose and write down five interesting candi- date models that you could build with these predictors.
Candidate predictors for the model are: gender: 0=male, 1=female age: age in years educ: educational level, 0=more than primary, 1=only primary consult: number of consults in the last 12 months household: 0=with others, 1=alone partner: 0=with partner, 1=single diagnosis: wether GP made diagnosis, 0=diagnosis , 1=only complaint (no diagnosis) problems: presentation of mental or social problems at inclusion, 0=no, 1=yes somatic: presentation of somatic complaints at inclusion, 0=no, 1=yes antidep: antidepressant prescribed in the last 12 months, 0=no, 1=yes depGP: depression code at GP, 0=no, 1=yes cidi1: any depressed feelings in lifetime, 0=no, 1=yes cidi2: any loss of interest in lifetime, 0=no, 1=yes
All could be important, I choose:
models:
First inspect the data
str(dat_prac2)
'data.frame': 1046 obs. of 14 variables:
$ depr : num 0 0 1 0 0 0 0 0 0 0 ...
$ gender : num 0 1 0 1 1 0 1 1 1 1 ...
$ age : num 66 56 51 24 60 49 54 57 63 52 ...
$ consult : num 2 15 24 7 6 12 17 10 8 1 ...
$ educ : num 0 0 1 0 0 0 0 0 0 0 ...
$ household: num 0 0 0 0 0 0 0 0 0 1 ...
$ partner : num 0 0 0 1 0 0 0 0 0 1 ...
$ diagnosis: num 0 0 0 1 1 0 1 0 1 1 ...
$ problems : num 0 0 0 0 1 0 0 0 0 0 ...
$ somatic : num 1 1 1 1 0 1 1 1 1 1 ...
$ depGP : num 0 0 0 0 0 0 0 0 0 0 ...
$ antidep : num 0 0 0 0 0 0 0 0 0 0 ...
$ cidi1 : num 0 0 0 0 0 1 1 1 0 0 ...
$ cidi2 : num 0 0 1 0 0 1 0 0 0 0 ...
Setup a container for the models
library(formula.tools)
models <- data.frame(
model = 0:4,
outcome = "depr",
rhs_char = c("1",
"gender",
"gender + age",
"gender + age + problems",
"gender + age + problems + household")
)
models %<>%
mutate(
model_file = paste0("bayes_2_Model2.", model, ".txt"),
model_file_full = here("analysis", model_file),
formula_char = paste0(outcome, "~", rhs_char),
covariates = map(formula_char, ~rhs.vars(formula(.x))),
parameters = map(covariates, ~{
params = c("alpha")
if (length(.x) > 0) params = c(params, paste0("b.", .x))
return(params)})
)
models$parameters
[[1]]
[1] "alpha"
[[2]]
[1] "alpha" "b.gender"
[[3]]
[1] "alpha" "b.gender" "b.age"
[[4]]
[1] "alpha" "b.gender" "b.age" "b.problems"
[[5]]
[1] "alpha" "b.gender" "b.age" "b.problems" "b.household"
Declare starting values for each model
These should be lists with values named per variable
We will set starting values between -5 and 5 for each of them, as they are all regression coefficients. To make sure not all variables start with the same coefficient, we will cycle the starting values
nmodels <- nrow(models)
start_values <- c(-5, 0, 5)
nchains = length(start_values)
# models$formula_char %>%
# map(formula) %>%
# map(rhs.vars) %>%
# map(print)
start_values_list <- vector(length = nmodels, mode = "list")
for (i in 1:nmodels) {
params = models$parameters[i][[1]]
nparams = length(params)
start_values_model <- vector(length = nchains, mode = "list")
for (j in 1:nchains) {
start_values_model[[j]] <- vector(length = nparams, mode = "list")
names(start_values_model[[j]]) <- params
## assign value of alpha
start_values_model[[j]][["alpha"]] <- start_values[j]
## assign values of other variables
if (nparams > 1) {
for (k in 2:(nparams)) {
shifter = 1 + ((1 + j + k) %% 3)
start_values_model[[j]][[k]] <- start_values[shifter]
}
}
}
start_values_list[[i]] <- start_values_model
}
print(start_values_list[[2]])
[[1]]
[[1]]$alpha
[1] -5
[[1]]$b.gender
[1] 0
[[2]]
[[2]]$alpha
[1] 0
[[2]]$b.gender
[1] 5
[[3]]
[[3]]$alpha
[1] 5
[[3]]$b.gender
[1] -5
Compile the models, store them in a list
start_values_list[[2]][[1]]
$alpha
[1] -5
$b.gender
[1] 0
dat_prac2_cage <- mutate(dat_prac2, age = (age - min(age)) / diff(range(age)))
j_models <- map2(models$model_file_full, start_values_list,
~jags.model(.x,
data = dat_prac2_cage,
inits = .y,
n.chains = 3))
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "gender" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "age" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "consult" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "educ" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "household" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "partner" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "diagnosis" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "problems" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "somatic" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "depGP" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "antidep" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi1" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 1
Total graph size: 1051
Initializing model
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "age" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "consult" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "educ" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "household" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "partner" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "diagnosis" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "problems" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "somatic" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "depGP" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "antidep" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi1" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 2
Total graph size: 2104
Initializing model
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "consult" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "educ" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "household" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "partner" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "diagnosis" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "problems" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "somatic" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "depGP" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "antidep" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi1" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 3
Total graph size: 3388
Initializing model
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "consult" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "educ" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "household" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "partner" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "diagnosis" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "somatic" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "depGP" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "antidep" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi1" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 4
Total graph size: 4566
Initializing model
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "consult" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "educ" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "partner" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "diagnosis" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "somatic" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "depGP" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "antidep" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi1" in data
Warning in jags.model(.x, data = dat_prac2_cage, inits = .y, n.chains = 3):
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 5
Total graph size: 5772
Initializing model
models %<>%
mutate(jags_model = j_models,
init_values = start_values_list,
out_vars = c("alpha", setdiff(rhs_char, "1")))
Do not forget to specify prior distributions and initial values for all parameters in each logistic regression model you run. For each model, inspect the convergence, examine and write down the posterior results: posterior means and 95% CIs of the parameters in each model, and the DIC.
Now burnin each model and extract the relevant parameters
nburnin = 2000
niter = 2000
burnin_and_run <- function(model, params, nburnin, niter) {
update(model, nburnin)
# out_vars <- c("alpha")
# if (length(vars) > 0) {
# out_vars <- c(out_vars, paste0("b.", vars), paste0("or.", vars))
# }
samples = coda.samples(model, n.iter = niter, variable.names = params)
return(samples)
}
models %<>%
mutate(samples = map2(jags_model, parameters,
~burnin_and_run(.x, .y, nburnin = nburnin, niter = niter)))
Check convergence for the models
(implement in ggplot later)
# models %>%
# mutate(sample_df = samples %>%
# map(~as.data.frame(.x[[1]])) %>%
# rbindlist(idcol = "chain", fill = T))
plot(models$samples[1][[1]])
plot(models$samples[2][[1]])
plot(models$samples[3][[1]])
plot(models$samples[4][[1]])
plot(models$samples[5][[1]])
Extract the DIC
models %<>%
mutate(dic = map(jags_model, ~dic.samples(.x, n.iter = niter, type = "pD")))
Look at DICs
get_dic <- function(x) {
deviance <- sum(x$deviance)
psum <- sum(x[[2]])
penalized_deviance = deviance + psum
return(penalized_deviance)
}
models %<>%
mutate(dic_value = map_dbl(dic, get_dic))
models %>%
select(outcome, rhs_char, dic_value)
outcome rhs_char dic_value
1 depr 1 886.6655
2 depr gender 881.7402
3 depr gender + age 881.6680
4 depr gender + age + problems 873.2346
5 depr gender + age + problems + household 870.8503
Obtaining the right variables for your analysis and putting these data in a list For each model you need to specify the right data. For this you need to see in which columns of the dataframe dat_prac2 your variables of interest are stored. You need to select these columns, store them in an object, and then store those in a list for rjags. For example, if you want to select the variables depr, gender and partner, you can go about it in the following way:
head(dat_prac2) ###Use this code to see an abbreviated version of the dataframe, and count in which columns your variables are ### we see that depr is in column 1, gender in 2, and partner in 7. ##Store each column you need in an object with an appropriate name: depr <- dat_prac2[,1] gend <- dat_prac2[,2] part <- dat_prac2[,7]
n=length(data_prac2_ex1[,1]) ### don’t forget to specify n too - the number of observations which is needed for your for loop.
Now state which object belongs with which variable in your jags model code, and put this in a list you can supply to rjags:
data_m1 <- list(“n”=n, “depr”=depr, “gender” = gend, “partner” = part) Monitoring the DIC To keep track of the dic, add the following (adjusted) line of code after using coda.samples()
in the first argument you specify which rjags model to sample the dic for, the second how many samples to take, and the third the type of penalty function (pD for the DIC).
dic.mod1 <- dic.samples(model = rjags_ex1, n.iter=2000, type = “pD”) To obtain the dic, call dic.mod1 and write down the penalized deviance.
All
Report and interpret the posterior results of your final model.
Exercise 3: Prediction with a Bayesian Logistic Regression Model In this exercise you will predict the probability of depression for a new patient using Bayesian prediction.
Add the data for the predictor variables for the new patient at the bottom (last row) of the data frame dat_prac2. Specify the depression score as NA, indicating that this is a missing observation.
You can do this with the following code:
load(“data_practical2.Rdata”)
data_new_patient=c(NA,2,3.50000E+01,2.00000E+00,0.00000E+00,0.00000E+00, 0.00000E+00,1.00000E+00,1.00000E+00,1.00000E+00,
0.00000E+00,0.00000E+00,1.00000E+00,1.00000E+00)
dat_prac2[1047,] <- data_new_patient
dat_prac2_cage_pred <- mutate(dat_prac2, age = (age - min(age)) / diff(range(age)))
Is the new patient male or female? Has he/she been prescribed antidepressants?
Report and interpret the 95% CI of the predicted probability of depression for this patient. Don’t forget to request depr[1047] and p[1047] as parameters to save samples for. What are your conclusions?
final_init <- start_values_list[[5]]
final_params <- models$parameters[5][[1]]
jm_predict <- jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"),
data = dat_prac2_cage_pred,
inits = final_init, n.chains = length(final_init))
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "consult" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "educ" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "partner" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "diagnosis" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "somatic" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "depGP" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "antidep" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "cidi1" in data
Warning in jags.model(here("analysis", "bayes_2_Model2.4_predict.txt"), :
Unused variable "cidi2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 1046
Unobserved stochastic nodes: 6
Total graph size: 5780
Initializing model
update(jm_predict, nburnin)
samp_predict <- coda.samples(jm_predict, variable.names = c(final_params, "pred_depr", "pred_p"),
n.iter = niter)
summary(samp_predict)
Iterations = 3001:5000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha -2.0474 0.2601 0.003358 0.016133
b.age -0.3861 0.3219 0.004155 0.016242
b.gender 0.4826 0.1974 0.002548 0.007641
b.household 0.4694 0.2202 0.002842 0.004432
b.problems 0.7252 0.2465 0.003182 0.004746
pred_depr 0.3823 0.4860 0.006274 0.006440
pred_p 0.3837 0.0828 0.001069 0.002342
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha -2.55852 -2.2240 -2.0386 -1.8701 -1.5495
b.age -1.01142 -0.6049 -0.3912 -0.1674 0.2481
b.gender 0.08824 0.3525 0.4801 0.6137 0.8647
b.household 0.03785 0.3239 0.4712 0.6120 0.9038
b.problems 0.22526 0.5631 0.7259 0.8932 1.1976
pred_depr 0.00000 0.0000 0.0000 1.0000 1.0000
pred_p 0.23241 0.3252 0.3813 0.4387 0.5536
Compare this Bayesian approach of handling missing data with classical missing data methods you are familiar with. What are differences and similarities?
Does not work with missing values
Today’s considerations for your presentation on Friday: Fit a Bayesian logistic regression model on your own data
Discuss differences and similarities between the DIC and other information criteria.
Use a Bayesian model for prediction with your own data set
Discuss missing data analysis from a Bayesian perspective
sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.6
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] bindrcpp_0.2 formula.tools_1.7.1 Bain_0.1.0
[4] rjags_4-6 coda_0.19-1 broom_0.4.2
[7] epistats_0.1.0 ggplot2_2.2.1 here_0.1
[10] purrr_0.2.4 magrittr_1.5 data.table_1.10.4
[13] dplyr_0.7.4
loaded via a namespace (and not attached):
[1] Rcpp_0.12.14 git2r_0.20.0 plyr_1.8.4
[4] bindr_0.1 class_7.3-14 tools_3.3.2
[7] digest_0.6.13 evaluate_0.10.1 tibble_1.3.4
[10] gtable_0.2.0 nlme_3.1-131 lattice_0.20-35
[13] pkgconfig_2.0.1 rlang_0.1.6 psych_1.7.5
[16] yaml_2.1.16 parallel_3.3.2 mvtnorm_1.0-6
[19] e1071_1.6-8 stringr_1.2.0 knitr_1.18
[22] rprojroot_1.2 grid_3.3.2 glue_1.2.0
[25] R6_2.2.2 foreign_0.8-69 rmarkdown_1.8
[28] tidyr_0.7.2 reshape2_1.4.2 operator.tools_1.6.3
[31] MASS_7.3-47 backports_1.1.0 scales_0.4.1
[34] htmltools_0.3.6 assertthat_0.2.0 mnormt_1.5-5
[37] fungible_1.5 colorspace_1.3-2 nleqslv_3.3.1
[40] labeling_0.3 stringi_1.1.6 lazyeval_0.2.0
[43] munsell_0.4.3
This R Markdown site was created with workflowr


Change the validity state, and notifies the form when the control changes validity. (i.e. it does not notify form if given validator is already marked as invalid).
This method should be called by validators - i.e. the parser or formatter functions.
validationErrorKey – {string} –
Name of the validator. the validationErrorKey
will assign
to $error[validationErrorKey]=isValid
so that it is available for data-binding.
The validationErrorKey
should be in camelCase and will get converted into dash-case
for class name. Example: myError
will result in ng-valid-my-error
and ng-invalid-my-error
class and can be bound to as {{someForm.someControl.$error.myError}}
.
isValid – {boolean} –
Whether the current state is valid (true) or invalid (false).
Read a value from view.
This method should be called from within a DOM event handler.
For example input
or
select
directives call it.
It internally calls all formatters
and if resulted value is valid, updates the model and
calls all registered change listeners.
value – {string} –
Value from the view.
Actual string value in the view.
The value in the model, that the control is bound to.
Whenever the control reads value from the DOM, it executes all of these functions to sanitize / convert the value as well as validate.
Whenever the model value changes, it executes all of these functions to convert the value as well as validate.
An bject hash with all errors as keys.
True if user has not interacted with the control yet.
True if user has already interacted with the control.
True if there is no error.
True if at least one error on the control.


[[macros.html]]
Register Correction Application
[[common/style.html]]
[[correction/businessDetails.html]]
[[change-of-registration/legal-name.html]]
[[change-of-registration/nature-of-business.html]]
[[change-of-registration/addresses.html]]
{% if business.legalType == 'GP' %}
[[change-of-registration/partner.html]]
{% elif business.legalType == 'SP' %}
[[change-of-registration/proprietor.html]]
{% endif %}
{% if correction.type == 'CLIENT' %}
[[change-of-registration/completingParty.html]]
{% endif %}


If you have found a vulnerabillity in lolbot, please send it to any of my contact methods (even if it's my email, or DMs, or IRC PMs.)
PLEASE DO NOT:
Report vuln in public with reproduction instructions
This will cause trouble users to exploit the bot and thus ruining the security of the software
- even if you didn't mean to.
Use it against the software
This will most likely result in something bad.
Thank you!
keyboard_arrow_left


Create a folder if necessary by writing the word math at the "folder to open" cue. When asked if you want to create the folder, press y for yes. Make a new file within the folder called "practice 1" for use during this lesson. Again, just write the file name when prompted to do so. Answer y for yes when asked if you want to create the file. Give your student help as needed.
in-context


New Orleans after Katrina. Photo courtesy Digital Globe. Panning/Zooming software courtesy of Michal Migurski and his GSV software. Thanks to Rafe Colburn for pointing out the picture. Among others, Mennonite Disaster Service is seeking donations for disaster response.


/customer/feedback/update
Update feedback record.
Parameters
| Name |
Type |
Required |
Example |
Description |
| clientId |
string(40) |
yes |
04f78...9f9f3 |
Client id is a unique user identifier provided by GetFiveStars. |
| customerCustomId |
string |
no |
123 |
Customer custom id. |
| customerId |
integer |
yes |
1 |
Customer id. |
| customerJobId |
string |
no |
123 |
Customer job id. |
| hash |
string(64) |
yes |
bd13a40...303f3c5 |
Hash is a digital request signature. Learn how to generate hash here. |
Request
POST https://getfivestars.com/api/customer/feedback/update
[Content-Type: application/json]
{
"customerId": 1,
"clientId": "04f78076fbfc66940169d2199af7c6698019f9f3",
"customerJobId": "1234",
"hash": "67fbca027648a3caf08e47ca685d772aa2d92fb526b8b6e91d74286301c99988",
}
Response
{
"errorCode": 0,
"errorMessage": "Success"
}
Errors
| errorCode |
errorMessage |
| -1 |
Unknown error |
| 1 |
Invalid hash sign |
| 2 |
Invalid clientId |
| 3 |
Server error - please contact support |
| 5 |
This user cannot manage selected business |
| 6 |
Invalid businessId |
| 44 |
Customer not found |
| 47 |
You can not update feedback history record |
| 48 |
Customer with this email already exists |


Umm, The Lenovo Smart Clock Is Cheaper Than Black Friday Right Now
Right now, the Lenovo Smart Clock is on sale for less than it will be on Black Friday, later this week. For Black Friday, the Lenovo Smart Clock is supposed to be $39.99 – which is half off of its regular price. Not bad, right? Well, Lenovo's official store on Rakuten has it for just $34.
To get this price, you'll need to sign in or sign up for a free Rakuten account (you get cash back, so it's definitely a good account to sign up for), and then use the promo code LEN6B at checkout.
The Lenovo Smart Clock is a very small smart display from the company, and it runs with the Google Assistant. It's meant to work as a smart display that sits next to your bed, showing you the time. But you can also ask Google Assistant for things like the weather and traffic, and of course, control your smart home products as well.
Advertisement
But it is sort of limited. You see, the Lenovo Smart Clock has a very small display, which means that it isn't all that great with playing video from YouTube, or even doing slideshows with your Google Photos albums. It's a smart display, but not as useful as one of the larger smart displays out there.
It is still really nice for putting on a bedside table, or even in your kids room. Since it is pretty small, and it does work well as a clock. But if you're getting this, planning to watch YouTube in bed on it, you'll be disappointed.
You can of course, use it to play music and such though, without any real issues.
Advertisement
You can pick up the Lenovo Smart Clock from Rakuten right now by clicking here, don't forget to use that code LEN6B at checkout to get it for just $34.


I was having a lot of trouble understand how JavaScript functions work and what everything means. The jargon was really confusing me and I was struggling to make progress. I ended up reading a whole bunch about it and watching YouTube videos to try get a better understanding. These things definitely helped but I felt like I needed to talk it over with someone and work through it at the same time. I ended up going and sitting with a great nerd friend who refuses to tell me answers flat out, but rather makes sure I understand the logic before we go through step by step. Sometimes I'll have a small epiphany while he's talking and end up the problem with my code or realise what I need to do and excitedly smash away at the keys before he can get to the grand finale. Rude, I know. I go through a crazy myriad of emotions when I get blocked.  and I just want to lie face-down on the floor for a few months to recuperate. Obviously, I get over these moments pretty quickly (sometimes it takes a little longer and some pretty intense meditating, though), and then I feel determined AF! From struggling with this stuff, I learned that I really do just need a break sometimes but the ideal length is nowhere near as long as my instincts tell me. If I force myself to try again after a few hours or the next day, I tend to get on a roll again pretty quick. Times that I've given myself a couple of days because it was stressing me out were totally unproductive and I just spent those days worrying. I need to work to get past the frustrated, wanting-to-give-up stage quicker and get back on the horse, so to speak.
Pseudocode - I'm haven't really formally done this yet but I think I have a fair understanding of it. It's kinda just the stuff I say out loud to myself. Trying something - I'm super confident doing this. I type a lot of nonsensical junk but now and then it works! Rubber ducky method - I haven't done this but I do talk to myself a lot... Reading error messages - I've been getting better with this and finding Developer Tools super handy! Console.logging - I have no idea how to use this. Definitely something I should look back on... Googling - My googling skills are pretty legit but sometimes they just don't cut it. I need to look at better ways to find stuff again. I always forget all the little tips and tricks! Asking your peers for help - I'm pretty confident doing this, especially considering how many people have offered to help! Asking coaches for help - I haven't really done this yet but would be happy to if I needed to. My friend has managed to absorb that so far. Improving your process with reflection - I haven't set formal reflection time but when I wake up in the morning I often do this and I end up feeling more motivated and sometimes having good ideas.


This page uses a ported JavaScript version of OpenBabel so that much more file format can be supported. Check it by click on the "Load Data" or "Save to File" button in viewer above.


Coming soon
HashResearch | We are a team of researcher, we build products that rule the future.
Contact us
- #57 & Near Nagarjuna Preuniversity, Yelahaka, Bangalore - 64
- +91<PHONE_NUMBER>
+91<PHONE_NUMBER> -<EMAIL_ADDRESS>Core Research
Featured On
© 2016 HashResearch Labs. All rights reserved.


You are about to buy:
The PayPal Wars for $65.00
Ship to:
5 Temasek Boulevard #09-01 Suntec Tower Five 038985 Singapore


|
SimpleITK
1.1.0
|
Segments structures in images based on a user supplied edge potential map. More...
#include <sitkGeodesicActiveContourLevelSetImageFilter.h>
Public Types | |
| typedef RealPixelIDTypeList | PixelIDTypeList |
| typedef GeodesicActiveContourLevelSetImageFilter | Self |
| Public Types inherited from itk::simple::ImageFilter< 2 > | |
| typedef ImageFilter | Self |
| Public Types inherited from itk::simple::ProcessObject | |
| typedef ProcessObject | Self |
Public Member Functions | |
| Image | Execute (const Image &image1, const Image &image2) |
| Image | Execute (const Image &image1, const Image &image2, double maximumRMSError, double propagationScaling, double curvatureScaling, double advectionScaling, uint32_t numberOfIterations, bool reverseExpansionDirection) |
| GeodesicActiveContourLevelSetImageFilter () | |
| double | GetAdvectionScaling () const |
| double | GetCurvatureScaling () const |
| uint32_t | GetElapsedIterations () const |
| Number of iterations run. More... | |
| double | GetMaximumRMSError () const |
| std::string | GetName () const |
| uint32_t | GetNumberOfIterations () const |
| double | GetPropagationScaling () const |
| bool | GetReverseExpansionDirection () const |
| double | GetRMSChange () const |
| The Root Mean Square of the levelset upon termination. More... | |
| Self & | ReverseExpansionDirectionOff () |
| Self & | ReverseExpansionDirectionOn () |
| Self & | SetAdvectionScaling (double AdvectionScaling) |
| Self & | SetCurvatureScaling (double CurvatureScaling) |
| Self & | SetMaximumRMSError (double MaximumRMSError) |
| Self & | SetNumberOfIterations (uint32_t NumberOfIterations) |
| Self & | SetPropagationScaling (double PropagationScaling) |
| Self & | SetReverseExpansionDirection (bool ReverseExpansionDirection) |
| std::string | ToString () const |
| virtual | ~GeodesicActiveContourLevelSetImageFilter () |
| Public Member Functions inherited from itk::simple::ImageFilter< 2 > | |
| ImageFilter () | |
| virtual | ~ImageFilter ()=0 |
| Public Member Functions inherited from itk::simple::ProcessObject | |
| virtual void | Abort () |
| virtual int | AddCommand (itk::simple::EventEnum event, itk::simple::Command &cmd) |
| Add a Command Object to observer the event. More... | |
| virtual float | GetProgress () const |
| An Active Measurement of the progress of execution. More... | |
| virtual bool | HasCommand (itk::simple::EventEnum event) const |
| Query of this object has any registered commands for event. More... | |
| ProcessObject () | |
| virtual void | RemoveAllCommands () |
| Remove all registered commands. More... | |
| virtual | ~ProcessObject () |
| virtual void | DebugOn () |
| virtual void | DebugOff () |
| virtual bool | GetDebug () const |
| virtual void | SetDebug (bool debugFlag) |
| virtual void | SetNumberOfThreads (unsigned int n) |
| virtual unsigned int | GetNumberOfThreads () const |
Private Types | |
| typedef Image(Self::* | MemberFunctionType) (const Image &image1, const Image &image2) |
Private Member Functions | |
| template<class TImageType > | |
| Image | ExecuteInternal (const Image &image1, const Image &image2) |
Private Attributes | |
| double | m_AdvectionScaling |
| double | m_CurvatureScaling |
| uint32_t | m_ElapsedIterations |
| double | m_MaximumRMSError |
| nsstd::auto_ptr< detail::MemberFunctionFactory< MemberFunctionType > > | m_MemberFactory |
| uint32_t | m_NumberOfIterations |
| double | m_PropagationScaling |
| bool | m_ReverseExpansionDirection |
| double | m_RMSChange |
Friends | |
| struct | detail::MemberFunctionAddressor< MemberFunctionType > |
Additional Inherited Members | |
| Static Public Member Functions inherited from itk::simple::ProcessObject | |
| static bool | GetGlobalDefaultDebug () |
| static void | GlobalDefaultDebugOff () |
| static void | GlobalDefaultDebugOn () |
| static void | SetGlobalDefaultDebug (bool debugFlag) |
| static void | GlobalWarningDisplayOn () |
| static void | GlobalWarningDisplayOff () |
| static void | SetGlobalWarningDisplay (bool flag) |
| static bool | GetGlobalWarningDisplay () |
| static void | SetGlobalDefaultNumberOfThreads (unsigned int n) |
| static unsigned int | GetGlobalDefaultNumberOfThreads () |
| static double | GetGlobalDefaultCoordinateTolerance () |
| Access the global tolerance to determine congruent spaces. More... | |
| static void | SetGlobalDefaultCoordinateTolerance (double) |
| Access the global tolerance to determine congruent spaces. More... | |
| static double | GetGlobalDefaultDirectionTolerance () |
| Access the global tolerance to determine congruent spaces. More... | |
| static void | SetGlobalDefaultDirectionTolerance (double) |
| Access the global tolerance to determine congruent spaces. More... | |
| Protected Member Functions inherited from itk::simple::ProcessObject | |
| virtual unsigned long | AddITKObserver (const itk::EventObject &, itk::Command *) |
| virtual itk::ProcessObject * | GetActiveProcess () |
| virtual void | OnActiveProcessDelete () |
| virtual void | onCommandDelete (const itk::simple::Command *cmd) SITK_NOEXCEPT |
| virtual void | PreUpdate (itk::ProcessObject *p) |
| virtual void | RemoveITKObserver (EventCommand &e) |
| Protected Member Functions inherited from itk::simple::NonCopyable | |
| NonCopyable () | |
| Static Protected Member Functions inherited from itk::simple::ImageFilter< 2 > | |
| static void | FixNonZeroIndex (TImageType *img) |
| Static Protected Member Functions inherited from itk::simple::ProcessObject | |
| template<class TImageType > | |
| static TImageType::ConstPointer | CastImageToITK (const Image &img) |
| template<class TPixelType , unsigned int VImageDimension, unsigned int VLength, template< typename, unsigned int > class TVector> | |
| static Image | CastITKToImage (itk::Image< TVector< TPixelType, VLength >, VImageDimension > *img) |
| template<class TImageType > | |
| static Image | CastITKToImage (TImageType *img) |
| static const itk::EventObject & | GetITKEventObject (EventEnum e) |
| template<typename T > | |
| static std::ostream & | ToStringHelper (std::ostream &os, const T &v) |
| static std::ostream & | ToStringHelper (std::ostream &os, const char &v) |
| static std::ostream & | ToStringHelper (std::ostream &os, const signed char &v) |
| static std::ostream & | ToStringHelper (std::ostream &os, const unsigned char &v) |
Segments structures in images based on a user supplied edge potential map.
where is image intensity and is the derivative of Gaussian operator.
This implementation allows the user to set the weights between the propagation, advection and curvature term using methods SetPropagationScaling() , SetAdvectionScaling() , SetCurvatureScaling() . In general, the larger the CurvatureScaling, the smoother the resulting contour. To follow the implementation in Caselles et al paper, set the PropagationScaling to (the inflation or ballon force) and AdvectionScaling and CurvatureScaling both to 1.0.
Definition at line 82 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Setup for member function dispatching
Definition at line 183 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
Define the pixels types supported by this filter
Definition at line 94 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
| typedef GeodesicActiveContourLevelSetImageFilter itk::simple::GeodesicActiveContourLevelSetImageFilter::Self |
Definition at line 84 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
virtual |
Destructor
| itk::simple::GeodesicActiveContourLevelSetImageFilter::GeodesicActiveContourLevelSetImageFilter | ( | ) |
Default Constructor that takes no arguments and initializes default parameters
| Image itk::simple::GeodesicActiveContourLevelSetImageFilter::Execute | ( | const Image & | image1, |
| const Image & | image2 | ||
| ) |
Execute the filter on the input images
| Image itk::simple::GeodesicActiveContourLevelSetImageFilter::Execute | ( | const Image & | image1, |
| const Image & | image2, | ||
| double | maximumRMSError, | ||
| double | propagationScaling, | ||
| double | curvatureScaling, | ||
| double | advectionScaling, | ||
| uint32_t | numberOfIterations, | ||
| bool | reverseExpansionDirection | ||
| ) |
Execute the filter on the input images with the given parameters
|
private |
|
inline |
Definition at line 128 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 120 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Number of iterations run.
This is a measurement. Its value is updated in the Execute methods, so the value will only be valid after an execution.
Definition at line 155 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 104 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inlinevirtual |
Name of this class
Implements itk::simple::ProcessObject.
Definition at line 166 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 136 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 112 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 148 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
The Root Mean Square of the levelset upon termination.
This is a measurement. Its value is updated in the Execute methods, so the value will only be valid after an execution.
Definition at line 163 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 144 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Set the value of ReverseExpansionDirection to true or false respectfully.
Definition at line 143 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 124 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 116 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 100 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 132 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 108 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
inline |
Definition at line 140 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
virtual |
Print ourselves out
Reimplemented from itk::simple::ProcessObject.
|
friend |
Definition at line 187 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 199 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 197 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 205 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 193 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 189 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 201 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 195 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 203 of file sitkGeodesicActiveContourLevelSetImageFilter.h.
|
private |
Definition at line 207 of file sitkGeodesicActiveContourLevelSetImageFilter.h.


Most people want to be good at the things they do. One of the most challenging aspects of this concept is figuring out how to obtain that “goodness” level. I want to become a better coder. I’ve taken a number of different approaches at reaching this elusive adjective, and this post is another crack at motivating myself to get there.
Here are a few of the things that I’ve done to try to get good at coding:
- Spent a month learning HTML, CSS, and Git through Treehouse videos.
- Enrolled in Portland Code School’s Javascript Immersion Course.
- Coded every day for 255 days and counting..
- Went to Meet-Ups and talked with people.
- Studied books (Javascript: The Good Parts, Bootstrap Book, Ducket Book)
- Made things: www.thejbac.com, www.IsabellaYoga.com. www.YoungCreativeSpace.com, www.getcityhacker.com, among others.
- Tracked my time, and committed to coding at least 20 hours per week, if not more.
- Did freelancing work for real clients.
These things have all helped me to get where I am in my coding journey, but recently, I realized that learning from other people in a structured classroom setting helps to motivate me, undertake projects that get me out of my comfort zone, and enables me to meet more cool people in the Portland tech community.
So, in May, I enrolled in Portland Code School’s Advanced Front-End Development Course.
This course entails two night classes per week over the course of four weeks. It’s structured so that we get out of it what we put into it, and due to that, we’ll be covering a lot of material.
Here is a brief outline of some of the things that we’ll be going over:
- Ajax
- Preprocessors like SASS and Stylus .
- JQuery
- JavaScript
- Going deeper into HTML and CSS
- More bootstrap!
- Less
- Package manager’s like Yeoman
- Gulp
- API's
- Browser-side storage
- Going Deeper into Chrome Developer Tools
- Selenium Browser testing
- Tuts Plus/Codecademy tutorials
I don’t think that I’ll be able to be incredible at all of these things over the course of a month, but I do know that I’ll be able to build a solid foundation over this time period. I can then use this to propel me to go deeper into the aforementioned tools, thereby continuing my coding journey. Thus, the plan going forward is to create Github Repos for becoming good at some of the above-mentioned tools, and to talk about them here, in a living, breathing blogpost. Here goes nuthin…
Project 1
Tools: JQuery, JS, Bootstrap
Repo: https://github.com/StephenJayYoung/CatWebsite
Description: This is a simple website that uses Bootstrap, JQuery, and some JS to display cats that are available for adoption.
Thoughts: Man, I forgot how to do some of the JS I learned in school! I needed to get back into that, and this project has been helpful for getting my feet wet again. I realized that I do have the skills to be good at this language, and working on this project has been very motivating in that sense. Working on it with a partner, and engagin in peer programming, has also been helopful, insightful, and refreshing!


TODO : put description here
Here are Direct Links to the Moodle Pages:
All assignments should be sent through the moodle. All information on credits and deadlines will be published there.| Autumn – 2018 | Autumn – 2017 |
|---|---|
|
| Swift | User Interface |
|---|---|
| Data Persistence | Concurrency&Networking |


About
Production
Regions
Distilleries
Scotch Now
Blah blah blah stuff about scotch and such.
Haven't found what your looking for?
Search


Broken Down Code with Comments
Below I walk you through the code, breaking it down into managable chuncks. This is a good way to learn what each bit of code actually does.
Tutorial Mode: With Comments - Without Comments
Jump to: HTML - CSS - JAVASCRIPT
The HTML
The first thing you want to do is create the html to contain the time calculator. When you pick a date, the time elapse or remaining will display within an empty paragraph tag. This tag can be placed anywhere on the page so make it fit your layout. Be sure to include an id to allow for formatting with CSS and data input through JS later on.
Next, you'll create the label and input for the date. Because different browsers display differently, it is advised to use a table calendar and text input rather than a date type input. Place this where ever you'd the input box to be.
Right after the input box you'll want to code your table calendar. You will see this table until you complete the CSS to hide it. Then you will use JavaScript to call it when the user clicks in the text box. The code looks a little confusing but is really simple once broken down.
Create a div with an id (I used "cal") to contain your table. This will be referenced later to hide the table until called. Inside the first div, create another div with id containing the “previous” text, allowing users to go back a month. Below the previous div, but inside the first, create another div with id containing the “next” text, allowing the calendar to move forward a month.
Create a table with id. Inside the table, create an empty caption with id. This, along with the rest of the table, will be populated with JS.
Inside the table, you’re going to create a cell row using the tr tag. Inside the this first row you’ll make 7 empty header cells using the th tag. Inside each of these put the day of the week, starting with Sun.
Underneath the first table row (under the closing tr tag), you'll create another table row. Make 7 empty table cells using the td tag. Repeat step 7 four more times. You should have 5 rows with 7 cells in each. This does not include the row with headers.
After the closing table tag, and before the last ending div tag, create a paragraph with id. Inside this paragraph write “close”.
Jump to: HTML - CSS - JAVASCRIPT
The CSS
And that's it for the HTML portion of this tutorial. Next, we'll be creating the CSS to style and hide the calendar until clicked. The code below targets the div with the id “cal” and tells it what size it should be, the amount of spacing it should have, and then centers it. It also hides the calendar with the display: none attribute.
This next portion sets the “previous” text on the calendar to the left. And the “next” text to the right.
This targets the table tag with the “calBody” id. First, we need to clear the two floats above or the calendar won’t look right. Then we set the height and width to take up the entire div it’s contained within. We also center the table using the margin attribute.
Almost finished with the CSS, this targets the closing text and sets it in the bottom right corner of the “cal” div.
Jump to: HTML - CSS - JAVASCRIPT
The JavaScript
Here is where things get interesting. Everything below here is what makes the calculator function. You could say this is the background code, the engine. We’re going to start by telling the browser to render the code in strict mode. Then we’ll declare our global variables. The variables are visible to all functions and can be called at any time. As you read, please remember the text in green is a comment I added to further explain certain lines.
Below that, we’re going to create a function to display and populate our calendar. Inside that function we will begin by creating our local variables. Local variables can only be used within the function holding them. The monthArray variable will be cycled through via the previous and next text button in the tables HTML. The thing to remember about arrays is the first option is always 0. In this case January would be zero while march would be two.
Now that the local variables are done, we need to make the previous and next buttons actually work. When you click the previous button, this code tells the function to go back by one in the monthArray above. Otherwise, it tells it to move forward by one if the next button is clicked.
We’ve already listed our local variables, now it’s time to set their values. The code below pulls the data from the user’s selections within the table. The last tine in this section will write the captionValue into the empty captions we wrote above.
This next section is kind of large but it’s purpose is simple. It sets the number of days in each month to the appropriate count. All the months with 31 days are listed first by their number in the array. Remember that January is 0. The else if section sets up February with it’s unique numbers. If a leap year, it gets 29, if not it gets 28. Then follows the rest of the months set to 30 days.
Our next task is to populate the table cells with the calendar dates. We are also going to clear the cells of any current content. This allows new numbers to be inserted when the month is changed.
We’re going to need to count through the cells as we go so we set the current number to the “I” variable. This allows the value to change as the function continues looping. Aside from populating the cells with the day number, it also applies either the “date” class or the “futuredate” class which you can chose to style in the CSS. I chose not to style them for this example. The final task in this segment is the display the calendar. It changes the CSS to display: block instead of none.
Here we will begin a new function to allow the user to select a date by clicking in the table and to store their selection. First we need to have the function execute when a page event happens (when clicked).
This next chunk keeps the calendar open if nothing is selected and also keeps the input box empty. Now we need to gather the user’s current date. Having their current date allows the code to calculate has long has passed or how long until the date chosen by the user. The last thing to do in this function is to close the calendar table and to call the countdown function.
Once the function has done its job and the user has chosen a date, it’s time to close the calendar and update the countdown. This will be done with two separate functions, so we’ll add some code to run those functions at the end of the current one.
Next, we need to create the function to hide the calendar. This is a simple one that changed the CSS back to none.
These two functions below cycle backwards and forwards through the months as the user clicks previous and next.
Now it’s time to use the date selected above to populate the countdown. Our first job is to create the local variables to store the current date and the date chosen by the user. These will be used by the “if” and “else if” statements in the next two segments.
Once the functions variables are defined, we need to create the first statement which will calculate the time until a future date. To do this we will start with the years, so we’re going to subtracts the future date from the current one, and then divide that number by the number of milliseconds in a year. You can see this behind done in the variable “yearsUntil.” Once the math is done, we need to write the results to the paragraph with the id “elapTime.”
Note that any text you want written needs to be placed in quotes. To place a variable value inside direct text, you can concatenate by adding a + sign as shown below. Any spaces need to be placed inside the quotes as well.
We’re going to repeat the step above for the months and days, but this time we’re going to divide the years until by the number of milliseconds in a month and then in a day.
Directly after the above “if” statement, we need to add an “if else” to calculate the time passed form the time chosen. This will look almost identical to the last statement except instead of subtracting the chosen date from today, we’re going to subtract today from the chosen date. This will give us the time passed from a past date.
In order for the functions above to work, they need to be called when a chosen event happens. Such as clicking a button. To make this happen we will create the event listeners to fire when their selected elements are clicked on. This is a three-step process for each desired event. First is a variable to identify the element we want the function to “listen” to. These elements are selected either by the assigned id or HTML tag type. The we create an “if” statement to fire when it’s clicked. This will cause it to call the desired function. This is followed by an “else if” statement that does the same thing for older browsers.
This next event listener is similar to the first, however it needs to be able to tell which cell in the calendar is clicked. To allow this we need to create a for statement to count the cells when clicked. Remember to add this for statement to both the if and the else if statements.
There are still two more event listeners needed. One to listen for when the close button is clicked on one to listen to the previous and next buttons. Because the second is listening for two events, it will need two variables and instructions for both.
Our final step is to create a listener to call the event listener function above. This one will execute as soon as the page finishes loading.


Interview Listings
Display listings by programming language:
None
Java
Python
C++
C#
C
JavaScript
Refresh


If you want to write your card up by hand, you can find the syntax & structure below.
But for the sake of convenience you can use the following generator.
Copy the code above and paste it into a new gist as described here.
Additionally, you can preview your card.
For the techy people!
| Line no. | Field | Example |
|---|---|---|
| 1 | Name | Jonas Kohl |
| 2 | Pronoun(s) (multiple seperated by or) |
he/him or they/them |
| 3 | Bio paragraph | Hello, this is me! |
| 4 | Links starting delimiter (has to be links: ) |
links: |
| 5+ | Links (in the format site@ ) |
github@jonaskohl |
| …+1 | Links ending delimiter (has to be end; ) |
end; |
| Every field below is optional1 | ||
| …+2 | Foreground color (valid CSS color) | #000000 |
| …+3 | Background color (valid CSS color) | #ffffff |
| …+4 | Link color (valid CSS color) | #0000ff |
| …+5 | Link color (hovered & focused) (valid CSS color) | #ff0000 |
1If a field is defined, every field before also needs to be defined (for example, you can only set the link color if both the foreground and background color are defined)


| name_spec {vctrs} | R Documentation |
A name specification describes how to combine an inner and outer names. This sort of name combination arises when concatenating vectors or flattening lists. There are two possible cases:
Named vector:
vec_c(outer = c(inner1 = 1, inner2 = 2))
Unnamed vector:
vec_c(outer = 1:2)
In r-lib and tidyverse packages, these cases are errors by default, because there's no behaviour that works well for every case. Instead, you can provide a name specification that describes how to combine the inner and outer names of inputs. Name specifications can refer to:
outer
: The external name recycled to the size of the input
vector.
inner
: Either the names of the input vector, or a sequence of
integer from 1 to the size of the vector if it is unnamed.
name_spec, .name_spec |
A name specification for combining
inner and outer names. This is relevant for inputs passed with a
name, when these inputs are themselves named, like
See the name specification topic. |
# By default, named inputs must be length 1: vec_c(name = 1) # ok try(vec_c(name = 1:3)) # bad # They also can't have internal names, even if scalar: try(vec_c(name = c(internal = 1))) # bad # Pass a name specification to work around this. A specification # can be a glue string referring to `outer` and `inner`: vec_c(name = 1:3, other = 4:5, .name_spec = "{outer}") vec_c(name = 1:3, other = 4:5, .name_spec = "{outer}_{inner}") # They can also be functions: my_spec <- function(outer, inner) paste(outer, inner, sep = "_") vec_c(name = 1:3, other = 4:5, .name_spec = my_spec) # Or purrr-style formulas for anonymous functions: vec_c(name = 1:3, other = 4:5, .name_spec = ~ paste0(.x, .y))


Installation and Setup¶
Dependencies¶
Aqua Chemistry is built upon Aqua. Like Aqua, at least Python 3.5 or later is needed to use Qiskit Aqua Chemistry. In addition, Jupyter Notebook is recommended for interacting with the tutorials. For this reason we recommend installing the Anaconda 3 Python distribution, as it comes with all of these dependencies pre-installed.
Code Installation¶
We encourage you to install Aqua Chemistry via the pip package management system:
pip install qiskit-aqua-chemistry
pip will handle all dependencies automatically (including the dependencies on Aqua and Qiskit Core). and you will always install the latest (and well-tested) release version.
If your intention is not so much to access Aqua Chemistry as a tool to perform chemistry computations on a quantum machine, but rather to extend Aqua Chemistry with new research contributions — such as new algorithms, algorithm components, input-translation operators or drivers — then it is advisable to clone both the Aqua Chemistry and Aqua Git repositories in order to have easier access to the source code of the various components.
Note
We recommend using Python virtual environments to improve your experience.
Jupyter Notebooks and input files for Aqua Chemistry are included as part of the Aqua Tutorials.
Installation of Chemistry Drivers¶
To run chemistry experiments on various molecules, you will also need to install one of the supported classical computational chemistry programs, or drivers, interfaced by Aqua Chemistry. Currently, Aqua Chemistry comes with built-in interfaces for four drivers:
- Gaussian™ 16, a commercial chemistry program
- PSI4, an open-source chemistry program built on Python
- PySCF, an open-source Python chemistry program
- PyQuante, a pure cross-platform open-source Python chemistry program
While the logic to
interface these drivers is supplied as part of the Aqua Chemistry installation, the dependent chemistry programs
need to be installed separately. This can be done by following the instructions provided in Section “Drivers”.
Supporting additional drivers in Aqua Chemistry can be easily achieved by extending the BaseDriver
interface.
Even without installing any of the drivers above, it is still possible to run chemistry experiments by passing to the inout-translation layer a Hierarchical Data Format 5 (HDF5) binary file serializing the intermediate data previously generated by one of the supported chemistry drivers. This offers researchers the opportunity to share chemistry input files and replicate each other’s results. Given its support to take an HDF5 files as the input to initiate a chemistry experiment, AQUA Chemistry lists HDF5 as an additional driver — in fact, the only built-in driver coming with Aqua Chemistry.
A few sample HDF5 files are provided as input files in the chemistry
folder of the
Aqua Tutorials repository.


Yes, link-only answers are bad. I knew this even at the time
that I posted this answer, but I didn't feel like copying and
pasting was OK due to copyright, plagiarism, and what have you.
However, I now feel like it's OK since I have linked to the
original article. I should also mention the author's name,
though, for credit: Nicolas Gallagher. Here is the meat of the
article (note that "Thierry’s method" is referring to
----AFTER----


Issue 109
Welcome to another issue of Haskell Weekly!  This is a weekly summary of what’s going on in its community.
Want to see something featured in Haskell Weekly? Open an issue or pull request on GitHub.
Featured
-
The GHC team is pleased to announce the availability of GHC 8.4.3. Thanks to everyone who has contributed to developing, documenting, and testing this release!
-
Fibonacci compiles end-to-end: Haskell to WebAssembly via GHC
The project recently reached an exciting milestone: for the first time, simple Haskell functions can be compiled all the way down to WebAssembly from Haskell and executed with a real JavaScript runtime!
-
Combinatorics of permutations: Introduction
This is a very small chapter (just one page) in the book, and all we do here is define a permutation, introduce some notation, and count the total number of permutations of a given number of objects (everyone knows the answer).
-
Edit: A small package for rewriting things
I was working on a compiler in the last semester with some friends and we found it helpful in saving the “was something changed” information when doing constant folding (expression level) and passing that information to the Hoopl library.
-
Haskell/GHC symbol search cheatsheet
Several features of Haskell/GHC have low Google-ability because some of them are composed of symbols. This page is a reference collection to support search of them.
-
Even if you’re not interested in the music part, it’s an example of a Haskell program with a GUI, REPL based interface, tests, profiling, and a Shake-based build system, although probably idiosyncratic implementations of all of the above.
-
I’ve worked on a lot of open source Haskell projects over the past ten years, and continue to do so today. However, I’m making a conscious decision to shift where my focus lies. My goal is to focus much more on enabling others to work on projects, rather than doing so myself.
-
Squashed c x
let a library writer providex
in “c
-irrelevant” way to a library user. TheMonad
instance allows to work on the wrapped value. However, we cannot extract the original value. This restriction maybe be useful to enforce correctness, without relying on the module system! -
Looking back, the experience of upgrading Shake didn’t turn out to be as difficult as I had once thought. This exercise pretty much reinforced what I had already known about using Haskell industrially.
Jobs
-
Software engineer at ITProTV in Gainesville, Florida
ITPro.TV is a fast-growing digital media business that focuses on continuing education in technical domains. We are currently accepting applications for full-stack software professionals to join our small but talented multidisciplinary team.
In brief
- animate-frames: A workflow tool for converting sprite frames into a spritesheet
- Duet: A programming language focused on interactive collaboration between the developer and the computer
- Freq: A way to train a model that predicts the “randomness” of an input
- Primitive <IP_ADDRESS> released
- Refined: Refinement types with static and runtime checking
Package of the week
This week’s package of the week is sixty-five-oh-two, an eDSL for writing 65(C)02 bytecode.
Call for participation
Events
- May 31: GHC Newcomers Guide in Zürich, Switzerland
- June 2: Hang out, chat FP, work on some code in Boston, Massachusetts, United States
- June 3: Why Scala? in Naperville, Illinois, United States
- June 4: Haskell Peer Study Group in Vancouver, British Columbia, Canada
- June 5: Monadic Warsaw #14 in Warsaw, Poland
- June 6: Auckland FP Meetup: 3 Talks: Remote Haskell; Prototyping for Radar; Pure TensorFlow Combinators in Auckland, New Zealand
- June 7: Haskell.SG June Haskell Meetup in Signapore