model_samplesize.utf8

library(here)
library(MASS)
library(tidyverse)
require(rjags)
library(ggplot2)

This notebook uses JAGS to compute model predictions for the sampling frames experiment.

Prior distribution

First we set up the prior distribution over category means. We’re using a Gaussian process prior that favours smooth distributions.

# consider 7 spheres of different sizes
ncat <- 7
# suppose that the 7 spheres are equally spaced along the size dimension
test <- 1:ncat

# parameter specifying mean of the Gaussian process 
m <- 0
# parameters specifying covariance of the Gaussian process
sigma <- .5
tau <- 1.5
rho <- .1

mean_gp <- rep(NA, ncat)
cov_gp <- array(NA, c(ncat, ncat))

# mean and covariance matrix defining the Gaussian process
for(i in 1:ncat) {
  mean_gp[i] <- m
  cov_gp[i,i] <- (sigma^2) + (tau^2)
  if (i<7) {
    for(j in (i+1):ncat) {
      cov_gp[i,j] <- (tau^2) * exp(-rho * (test[i] - test[j])^2)
      cov_gp[j,i] <- cov_gp[i,j]
    }
  }
}

Let’s sample 6 functions from the prior to get a sense of what they look like.

# sample functions from the Gaussian process
nsample <- 6
samples <-  mvrnorm(n=nsample,mean_gp, cov_gp)

# define logistic transformation which maps real numbers to probabilities
logistic <- function(x) {
 return(1/(1+exp(-x)))
}

# apply the logistic transformation to the samples
category_means <- apply(samples, 1:2, logistic)

# plot the samples

category_means <- as_tibble(category_means) %>%
    set_names(c("1", "2", "3", "4", "5", "6", "7"))

## Warning: `as_tibble.matrix()` requires a matrix with column names or a `.name_repair` argument. Using compatibility `.name_repair`.
## This warning is displayed once per session.

category_means$samplenum <- paste("sample", 1:nsample)

tidymeans <- gather(category_means, cat, value, -samplenum)

pic <- tidymeans %>%
  ggplot(aes(x = cat, y = value, group=1)) +
  geom_line() +
  geom_point() +
  facet_wrap(~samplenum) +
  scale_y_continuous(lim=c(0,1)) +
  xlab("sizes") +
  ylab("prob of plaxium coating")
    
plot(pic)

Exercise

Run the chunk above a few times to draw samples from the prior. Does this prior capture the intuitive expectations that you bring to the spheres of Sodor problem? For example, the prior assigns some probability to U-shaped curves and other curves with one or more turning points: is this reasonable?

Sampling frames model

We’ll now use JAGS to model the generalization problem given to participants. There are two separate .bug files — one for the category sampling condition (category.bug) and one for the property sampling conditions (property.bug). The models for these conditions use the same prior but different likelihood functions, and the two .bug files reflect the two different assumptions about the likelihood.

# Script for running the sampling frames model using JAGS
# The code will give you some warnings about unused variables when you run it -- don't worry about them

make <- function(bugfile, obs = list()) {
  
  # parameters
  settings <- function(obs) {
    
    # by default there are 10 categories, although the experiment only
    # asks about 7
    if(!exists("ncat",obs)) obs$ncat <- 10
    
    # add locations for the categories
    obs$test <- 1:obs$ncat
    
    # add the dummy variables for the positive observations
    # (plaxium = 1)
    obs$plaxium <- rep.int(1,obs$nobs)
    
    # by default the prior over base rates is symmetric 
    # dirichlet with concentration .35
    if(is.null(obs$alpha)) {
      obs$alpha <- rep.int(.35,obs$ncat)
    }
    
    # default parameters for the gaussian process
    if(!exists("sigma",obs)) obs$sigma <- .5
    if(!exists("tau",obs)) obs$tau <- 1.5
    if(!exists("rho",obs)) obs$rho <- .1
    # curves will have a mean of 0.5
    if(!exists("m",obs)) obs$m <- 0 
     
    return(obs)
  }
  
  # initialise
  model <- list()
  
  # simulation parameters
  model$opt <- list(
    burnin = 20000,
    its = 100000,
    nchains = 1,
    thin = 10
  )

  # data to be given to JAGS
  model$obs = settings(obs)
  
  # store the jags model specification as a string
  model$string <- paste0(
    readLines(bugfile), 
    collapse="\n"
  )
  
  # construct the jags model object
  model$jagsmod <- jags.model(
    file = textConnection(model$string),
    n.adapt = model$opt$burnin,
    n.chains = model$opt$nchains,
    data = model$obs
  )
  
  # draw samples
  model$samples <- jags.samples(
    model = model$jagsmod, 
    variable.names = c("category_means"), 
    n.iter = model$opt$its,
    thin = model$opt$thin
  )
  
  # add a convenient summary
  model$out <- data.frame(
    test = model$obs$test,
    category_means  = apply(model$samples$category_means, 1, mean)
  )
  
  return(model)
}

sim <- list()

# --- simulations for sample size experiment---

sim$category_n2 <- make(
  bugfile = here("models","category.bug"),
  obs = list(nobs = 2, category = c(1,2))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "alpha" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 1
##    Total graph size: 161
## 
## Initializing model

sim$category_n6 <- make(
  bugfile = here("models","category.bug"),
  obs = list(nobs = 6, category = rep.int(c(1,2),3))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "alpha" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 6
##    Unobserved stochastic nodes: 1
##    Total graph size: 169
## 
## Initializing model

sim$category_n12 <- make(
  bugfile = here("models","category.bug"),
  obs = list(nobs = 12, category = rep.int(c(1,2),6))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "alpha" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 12
##    Unobserved stochastic nodes: 1
##    Total graph size: 181
## 
## Initializing model

sim$property_n2 <- make(
  bugfile = here("models","property.bug"),
  obs = list(nobs = 2, category = c(1,2))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "plaxium" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 2
##    Total graph size: 205
## 
## Initializing model

sim$property_n6 <- make(
  bugfile = here("models","property.bug"),
  obs = list(nobs = 6, category = rep.int(c(1,2),3))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "plaxium" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 6
##    Unobserved stochastic nodes: 2
##    Total graph size: 209
## 
## Initializing model

sim$property_n12 <- make(
  bugfile = here("models","property.bug"),
  obs = list(nobs = 12, category = rep.int(c(1,2),6))
)

## Warning in jags.model(file = textConnection(model$string), n.adapt =
## model$opt$burnin, : Unused variable "plaxium" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 12
##    Unobserved stochastic nodes: 2
##    Total graph size: 215
## 
## Initializing model

Plot the generalizations according to the models.

get_out <- function(mod, cond, samp) {
  out <- mod$out
  names(out)[2] <- "prediction" 
  out$sample_size<- samp
  out$condition <- cond
  return(out)
}

output <- rbind(
  get_out(sim$category_n2, "category", "small"),
  get_out(sim$property_n2, "property", "small"),  
  get_out(sim$category_n6, "category", "medium"),
  get_out(sim$property_n6, "property", "medium"),  
  get_out(sim$category_n12, "category", "large"),
  get_out(sim$property_n12, "property", "large") 
)

output<- output%>%
  mutate(sample_size = factor(sample_size, levels = c("small","medium","large"))) %>%
  # plot results only for the first seven categories along the size dimension
  filter(test <= 7)

pic2 <- output  %>%
  ggplot(aes(x = test, y = prediction, colour = condition)) +
  geom_line() +
  geom_point() +
  facet_wrap(~sample_size)
plot(pic2)

Day 4: Modeling category and property sampling

Day 4: Modeling category and property sampling

Charles Kemp

CHDSS 2019

Prior distribution

Exercise

Sampling frames model