Data Source

This analysis uses data in the paper by Milner et al on suicide by health professionals 2001 to 2012.

Data

Suicide by Profession

d1<-data.frame(read.csv("e:/pandoc/rmd/suicide01.csv",sep=","))
d2<-data.frame(read.csv("e:/pandoc/rmd/suicide02.csv",sep=","))

Data Manipulation:

library(dplyr)
fn<-function(occ){
  switch(as.character(occ),"Medical"=1,"Midwives and nurses"=2,"Other health professionals"=3,"Other occupations"=4)}
fg<-function(gender){switch(as.character(gender),"M"=1,"F"=2)}
d1$occupation<-lapply(d1$Occupation,fn)
d1$occupation<-as.integer(d1$occupation)
d1$gender<-lapply(d1$Gender,fg)
d1$gender<-as.integer(d1$gender)
d1<-d1%>%mutate(log_pop=log(Population))

Models

The following multilevel model incorporates the profession and gender as predictors for the suicide rate. The events rate is modeled as a Poisson process with log link function. Sampling was performed using rstan (Hamiltonian MCMC) as implemented in the rethinking package.

library(rethinking)
m.prof<-map2stan(
  alist(
    Suicides~dpois(lambda),
    log(lambda)<-a+profession[occupation]+sex[gender]+log_pop,
    a~dnorm(0,1),
    profession[occupation]~dnorm(0,1),
    sex[gender]~dnorm(0,1)
  ),
  data=d1,
  iter=1000
)
## 
## SAMPLING FOR MODEL 'Suicides ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.347 seconds (Warm-up)
##                0.376 seconds (Sampling)
##                0.723 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Suicides ~ dpois(lambda)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 0 seconds (Warm-up)
##                0 seconds (Sampling)
##                0 seconds (Total)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]
smp<-extract.samples(m.prof)
reg<-exp(smp$a+smp$profession[,1]+smp$sex[,1])/exp(smp$a+smp$profession[,1]+smp$sex[,2])
hist(reg,main="Relative suicide rate of medical practitioners: male to female",xlab="Male to female suicide rate")

p1m<-exp(smp$a+smp$profession[,1]+smp$sex[,1])/exp(smp$a+smp$profession[,2]+smp$sex[,1])
p2m<-exp(smp$a+smp$profession[,1]+smp$sex[,1])/exp(smp$a+smp$profession[,3]+smp$sex[,1])
p3m<-exp(smp$a+smp$profession[,1]+smp$sex[,1])/exp(smp$a+smp$profession[,4]+smp$sex[,1])
hist(p1m,main="Relative suicide rate in males",xlab="Medical to Midwives and nursing")

hist(p2m,main="Relative suicide rate: male ",xlab="Medical to Other health professionals")

hist(p3m,main="Relative suicide rate: male",xlab="Medical to Other occupations")

p1f<-exp(smp$a+smp$profession[,1]+smp$sex[,2])/exp(smp$a+smp$profession[,2]+smp$sex[,2])
p2f<-exp(smp$a+smp$profession[,1]+smp$sex[,2])/exp(smp$a+smp$profession[,3]+smp$sex[,2])
p3f<-exp(smp$a+smp$profession[,1]+smp$sex[,2])/exp(smp$a+smp$profession[,4]+smp$sex[,2])
hist(p1f,main="Relative suicide rate: female",xlab="Medical to Midwives and nurses")

hist(p2f,main="Relative suicide rate: female",xlab="Medical to Other health professionals")

hist(p3f,main="Relative suicide rate: female",xlab="Medical to Other occupations")

### Age and Suicicde Rates

Data preparation

af<-function(age){
  switch(as.character(age),
         "20-39"=1,
         "40-49"=2,
         "50-59"=3,
         "60-70"=4
         )
}
d2$age<-lapply(d2$Age,af)
d2$age<-as.integer(d2$age)
d2$gender<-lapply(d2$Gender,fg)
d2$gender<-as.integer(d2$gender)
d2$log_pop<-log(d2$Population)
m.age<-map2stan(
  alist(
    Suicides~dpois(lambda),
    log(lambda)<-a+sex[gender]+age_cat[age]+log_pop,
    a~dnorm(0,1),
    sex[gender]~dnorm(0,1),
    age_cat[age]~dnorm(0,1)
  ),
  data=d2,
  iter=1000
  
)
## 
## SAMPLING FOR MODEL 'Suicides ~ dpois(lambda)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)
##  Elapsed Time: 0.324 seconds (Warm-up)
##                0.282 seconds (Sampling)
##                0.606 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'Suicides ~ dpois(lambda)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 0 seconds (Warm-up)
##                0 seconds (Sampling)
##                0 seconds (Total)
## 
## [ 50 / 500 ]
[ 100 / 500 ]
[ 150 / 500 ]
[ 200 / 500 ]
[ 250 / 500 ]
[ 300 / 500 ]
[ 350 / 500 ]
[ 400 / 500 ]
[ 450 / 500 ]
[ 500 / 500 ]
precis(m.age,depth=2)
##             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a          -4.05   0.62      -4.94      -3.06    53 1.01
## sex[1]     -1.25   0.56      -2.13      -0.39    70 1.00
## sex[2]     -2.82   0.56      -3.67      -1.95    71 1.00
## age_cat[1] -0.94   0.46      -1.73      -0.32    62 1.00
## age_cat[2] -0.90   0.46      -1.69      -0.28    62 1.00
## age_cat[3] -1.15   0.46      -1.92      -0.52    63 1.00
## age_cat[4] -1.26   0.46      -2.11      -0.69    46 1.00
plot(precis(m.age,depth = 2))

sa<-extract.samples(m.age)
fun<-function(age,sex){
  exp(sa$a+sa$age_cat[,age]+sa$sex[,sex])
}
nam<-c("20-39","40-49","50-59","60-70")
boxplot(fun(1,1),fun(2,1),fun(3,1),fun(4,1),names=nam,xlab="Age",ylab="Suicide rate",main="Suicide rate by Age Category: Male")

boxplot(fun(1,2),fun(2,2),fun(3,2),fun(4,2),names=nam,xlab="Age",ylab="Suicide rate",main="Suicide rate by Age Category: Female")