Bayesian Modeling for human activity data (reported at the MathPsych2017)

author: Mr.Unadon@MrUnadon


Mac OS Sierra 10.12.1; R version3.3.2; rstan 2.10.1

Resut_ZIG_seed1234.png
Figure Zero-Inflated gamma dynamic model for human activity data.

1. For the atendee of MathPsych2017

Thank you for your visiting our page.

This article introduces tentative Bayesian models for human activity data.

If you have interest on the models introduced below, let’s discuss and develop the further suitable modeling.

Models has been proposed at the annual meeting of Mathematical Psychology 2017(Poster presentation).

I can send our poster to you (but only for those who attended the conference).

Author Imformation

logo.jpg
Name: Unadon(Certificate of licensed cook)
Skill: Bayesian Data Cooking

Interest and Future: Data Science for applied psychology and child welfare
Contact me : whitesnow1751★gmail.com
Twitter: @MrUnadon

2. Sample Data

2.1 Activity Data

I prepared a sample activity data (completly artificial and fictitious).

The one day activity data, N=3, time epocs = 144(recording activity per 10minutes in a day).

2.2 Read data

###################################################
# Bayesian Modeling for Human Activity Sample
###################################################

#clear work space
invisible({rm(list=ls());gc();gc()})

#library
library(tidyverse) # Here we use {dplyr} and {ggplot2}
library(rstan) #HMC
library(shinystan) #rstan results

#stan parallel
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())


#Read the completly artifical and fictitious sample data
Activation<-read.csv("https://mrunadon.github.io/images/SampleFictitiousActivation.txt",fileEncoding = "shift-jis")

#Stan Data
N=max(Activation$ID) #Sample 3ID 
t=c(ncol(Activation)-1)#144 time point
Act<-data.matrix(Activation[,-1])/100 #Activity Data Matrix
datastan = list(N=N,t=t,Act=Act)

Table Sample Data

2.3. Actograph

Visualizing the each person’s activity.

#data frame for ggplot2
dfActPlot<-Activation%>% #select base data
    tidyr::gather(key = TimePoint, value = Activity,-ID)%>% #reshape to long data
    dplyr::mutate(TimePoint=c(rep(1:(ncol(Activation)-1),each=nrow(Activation))), #integer time point
                  ID=as.factor(ID),
                  Activity=Activity/100) #factor ID

#Actograph
ggplot(dfActPlot,aes(x=TimePoint,y=Activity,colour=ID))+
    geom_line()+
    scale_x_continuous(expand = c(0,0))+
    facet_wrap(~ID,ncol = 1)

BDMactograph.png
Figure Actograph(raw data)

3. Gaussian Models

3.1. Overview

First, here I would like to introduce a popular state space model which is often called “Local Level Model”. This model estimates a successive activity state corresponding to a series of activity values. Moreover, it assumes the activity state moves with Gaussian noise epoch-by-epoch and observed values are generated from Gaussian distribution with that state as the expected value.

3.2. Model Code

R and Stan code are below.

#
# State: Gaussian, Observation: Gaussian Modeling
#
model_Gaussian<-"
data{
    int N;
    int t;
    matrix<lower=-1>[N,t] Act;
}


parameters{
    matrix <lower=0>[N,t]mu;
    real <lower=0>sigmaObs[N];
    real <lower=0>sigmaState[N];
    real <lower=0>muZero[N];
}

model{
    for (i in 1:N){
        sigmaObs[i] ~ cauchy(0,25);
        sigmaState[i] ~ cauchy(0,25);
        muZero[i] ~ normal(0,100);
        mu[i,1] ~ normal(muZero[i],sigmaState[i]);
        Act[i,1] ~ normal(mu[i,1],sigmaObs[i]);
        for (j in 2:t){
            mu[i,j] ~ normal(mu[i,(j-1)],sigmaState[i]);
            Act[i,j] ~ normal(mu[i,j],sigmaObs[i]);
        }
    }
}

"

#fit
fit = stan(model_code = model_Gaussian,
           data = datastan,
           chains = 4,
           iter = 4000,
           warmup = 2000,
           thin=2,
           seed=1234)

#summary function
mySummary=function(data){
    res<-data%>%
        data.frame()%>%
        tidyr::gather()%>%
        tidyr::separate(key,into = c("ID","TimePoint"), sep="\\.")%>%
        dplyr::mutate(ID=as.factor(stringr::str_replace(ID,pattern = "X",replacement = "")),
                      TimePoint=as.numeric(TimePoint))%>%
        dplyr::group_by(ID,TimePoint)%>%
        dplyr::summarize(EAP=mean(value),Upper95CI=quantile(value,0.975),Lower95CI=quantile(value,0.025))%>%
        dplyr::ungroup()%>%
        dplyr::left_join(dfActPlot,by=c("ID","TimePoint"))
    return(res)
}

#Apply Summarize Function
res=mySummary(rstan::extract(fit)$mu)


#Actograph
ggplot(res,aes(x=TimePoint,colour=ID,fill=ID))+
    geom_line(aes(y=Activity))+
    geom_line(aes(y=EAP))+
    geom_ribbon(aes(ymax=Upper95CI,ymin=Lower95CI),alpha=0.3)+
    facet_wrap(~ID,ncol = 1)+
    scale_x_continuous(expand = c(0,0))

3.3. Results(estimated state visualization)

Results showed the smooth state transitions. However, this model could not describe the “activity = zero”.

Result_Normal_seed1234.png
Figure Actograph and estimated state.

4. Dynamic Model with Exponential Distribution

4.1. Overview

Our second model also estimates a successive activity state corresponding to a series of activity values. However, it is assumed (1) the activity state moves with Gaussian noise, and (2) observed values are generated from exponential distribution with the activity state as the expected value.

4.2. Model Code

#
# State: Gaussian, Observation: Exponential Modeling
#
model_Exponential<-"
data{
    int t;
    int N;
    matrix<lower=0>[N,t] Act;
}


parameters{
    matrix <lower=0>[N,t]lambda;
    real <lower=0>sigmaState[N];
    real <lower=0>muZero[N];
}


transformed parameters {
    matrix <lower=0>[N,t]mu;
    for (i in 1:N){
        for (j in 1:t){
            mu[i,j] = 1/lambda[i,j];
        }    
    }
}
model{
    for(i in 1:N){
        lambda[i] ~ gamma(0.1,0.1);
        sigmaState[i] ~ cauchy(0,5);
        muZero[i] ~ normal(0,100);
        target+=normal_lpdf(mu[i,1]|muZero[i],sigmaState[i]);
        Act[i,1] ~ exponential(lambda[i,1]);
        for (j in 2:t){
            target+=normal_lpdf(mu[i,j]|mu[i,j-1],sigmaState[i]);
            Act[i,j] ~ exponential(lambda[i,j]);
        }
    }
}

"

#fit
fit2 = stan(model_code = model_Exponential,
           data = datastan,
           chains = 4,
           iter = 4000,
           warmup = 2000,
           thin=2,
           seed=1234)

#Apply Summarize Function
res2=mySummary(rstan::extract(fit2)$mu)

#posterior state
ggplot(res2,aes(x=TimePoint,colour=ID,fill=ID))+
    geom_line(aes(y=Activity))+
    geom_line(aes(y=EAP))+
    geom_ribbon(aes(ymax=Upper95CI,ymin=Lower95CI),alpha=0.3)+
    facet_wrap(~ID,ncol = 1)+
    scale_x_continuous(expand = c(0,0))

4.3. Results

Exponential distribution suitably described the zero value especially in sleeping time. However, this model indicated the MAP(maximum a positeriori) of activity observations are always zero. This assumption is not suitable for human activity.

Result_Exponential_seed1234.png
Figure Actograph and estimated state.

5. Dynamic Model for Binary data(1)

5.1. Overview

To estimate “expected sleep-awake probabilities”(we definitely need the validation with polysomnography or EEG. So here we call it “expected sleep-awake probability or state”), we composed the dynamic models for binary activity data. In the stan code, the activity data(originally non-negative) was transformed to binary data(zero or non-zero). The “expected sleep-awake state” is expressed by gaussian distribution. And estimated state are converted to bernoulli probability using inverse logit function.

5.2. Model Code

#
#
# Data: Binary transformed
# State: Gaussian, Observation: Bernoulli(inv_logit(Gaussian))
#
#
model_invLogit<-"
data{
    int N;
    int t;
    vector[t] Act[N]; //Activation record
}

parameters{
    matrix<lower=-5,upper=5> [N,t]mu;
    real <lower=0>sigmaState[N];
}

transformed parameters{
    matrix<lower=0,upper=1>[N,t]Mu;
    for(i in 1 : N){
        for(j in 1 : t){
            Mu[i,j] = inv_logit(mu[i,j]);
        }
    }
}

model{
    for (i in 1:N){
    int zero;
    int one;
        sigmaState[i] ~ normal(0,1);
        for (j in 2:t){
            mu[i,j] ~ normal(mu[i,(j-1)],sigmaState[i]);
            if (Act[i,j]==0){
                zero=0;
                target+=bernoulli_lpmf(zero|Mu[i,j]);
            }else{
                one=1;
                target+=bernoulli_lpmf(one|Mu[i,j]);
            }
        }
    }
}

"

datastan = list(N=N,t=t,Act=ifelse(Act==0,0,1))

#fit
fit4 = stan(model_code = model_invLogit,
            data = datastan,
            chains = 4,
            iter = 8000,
            warmup = 2000,
            thin=6,
            seed=1234)

#Apply Summarize Function
res4=mySummary(rstan::extract(fit4)$Mu)

#posterior state
ggplot(res4,aes(x=TimePoint,colour=ID,fill=ID))+
    geom_line(aes(y=EAP))+
    geom_point(aes(y=ifelse(res4$Activity==0,0,1)))+
    geom_ribbon(aes(ymax=Upper95CI,ymin=Lower95CI),alpha=0.3)+
    facet_wrap(~ID,ncol = 1)+
    scale_x_continuous(expand = c(0,0))

#for next modeling
pEAP=matrix(res4$EAP,ncol = 144,nrow=3,byrow = T)

5.3. Results

Result showed good transitions of “time series sleep-awake probability”.

These sucsessive probability will contribute to predict sleep or awake from actigraph.

Result_InvLogit_seed1234.png
Figure Binary transformed activity and the "expected sleep-awake probability".

6. Dynamic Model for Binary data(2)

6.1. Overview

Here I introduce another type of modeling, Bernoulli-Beta dynamic model. This model directly describe the “expected sleep-awake probability” using beta distribution. However, it is more difficult to interpret the estimated parameter values than introduced above. Moreover, this model will take a time about twice as large as avobe “inv_logit(gaussian) model”.

Now that we don’t recommend this modeling.

6.2. Model Code

#
# Data: Binary transformed
# (State): Beta, (Observation): Bernoulli
#

model_Beta<-"
data{
    int t; //data length per Sub
    int N; //N of Subject
    vector [t]Act[N]; //Activation record
}

parameters{
    vector <lower=0, upper=1> [t]thetaMu[N]; 
    real <lower=0.1> thetaPhi[N]; 
}
transformed parameters{
    vector <lower=0> [t] alpha[N];
    vector <lower=0> [t] beta [N];
    for(i in 1 : N){
        alpha[i] = thetaMu[i]*thetaPhi[i];
        beta [i] = (1-thetaMu[i])*thetaPhi[i];
    }
}
model {
    int zero;
    int one;
    zero=0;
    one=1;
    thetaPhi ~ gamma(0.01,0.01);
    thetaMu[1] ~ beta(1,1);
    for(i in  1: N){
        for(j in 2 : t ){
            thetaMu[i,j] ~ beta(alpha[i,j-1],beta[i,j-1]);
            if (Act[i,j]==0){
                target+=bernoulli_lpmf(zero|thetaMu[i,j]);
            }else{
                target+=bernoulli_lpmf(one|thetaMu[i,j]);
            }
        }
    }
}

"

datastan = list(N=N,t=t,Act=ifelse(Act==0,0,1))

#fit
fit3 = stan(model_code = model_Beta,
            data = datastan,
            chains = 4,
            iter = 8000,
            warmup = 2000,
            thin=6,
            seed=1234)

#Apply Summarize Function
res3=mySummary(rstan::extract(fit3)$thetaMu)

#posterior state
ggplot(res3,aes(x=TimePoint,colour=ID,fill=ID))+
    geom_line(aes(y=EAP))+
    geom_point(aes(y=ifelse(res3$Activity==0,0,1)))+
    geom_ribbon(aes(ymax=Upper95CI,ymin=Lower95CI),alpha=0.3)+
    facet_wrap(~ID,ncol = 1)+
    scale_x_continuous(expand = c(0,0))

6.3. Results

The “expected sleep-awake probability” is estimated smoothly with beta distribution. However, there is some theoretical unknowns (e.g. overdispersion and stochastic process) in this model.

Result_Beta_seed1234.png
Figure Binary transformed activity and the "expected sleep-awake probability".

7. Zero-Inflated Gamma modeling using estimated “sleep-awake probability”

7.1. Overview

This is our most chalenging model. Zero-Inflated Gamma dynamic modeling assumes (1) latent sleep-awake probability (already estimated with beta dynamic model), (2) activity state in both sleep and awake time, (3) activity fluctuations in both sleep and awake time, (4) state fluctuations in both sleep and awake time.

These latent parameters may contribute the developing for new insight in the psychological research using activity data.

If you have enough time epocs(minute-by-minute recording activity data=1440 time point per a day), sampling and thinning does not need so much like below.

7.2. Model Code

#
#
# Zero-Inflated Gamma Dynamic modeling
#
#

model_ZIG<-"
data{
    int t; 
    int N;
    vector [t]Act[N];
    vector [t]pEAP[N];
}

parameters{
    vector <lower=0,upper=10> [t] gammaMu [N];
    real <lower=0> gammaPhiSleep[N];
    real <lower=0> gammaPhiActive[N]; 
    real <lower=0> sigmaSleep[N];
    real <lower=0> sigmaActive[N];
}
transformed parameters{
    vector <lower=0> [t] shape[N]; 
    vector <lower=0> [t] rate [N]; 
    for(i in 1 : N){
        for(j in 1 : t){
            if (pEAP[i,j]<0.5){
                shape[i,j] = (gammaMu[i,j]^2)/gammaPhiSleep[i];
                rate[i,j]  = gammaMu[i,j]/gammaPhiSleep[i];
            }else{
                shape[i,j] = (gammaMu[i,j]^2)/gammaPhiActive[i];
                rate[i,j]  = gammaMu[i,j]/gammaPhiActive[i];
            }
        }
    }
}
model {
    real nearZero;
    nearZero=0.0000001;
    for(i in  1: N){
        gammaMu[i,1] ~ normal(0, 5);
        gammaPhiSleep[i]  ~ gamma(0.1, 0.1);
        gammaPhiActive[i]  ~ gamma(0.1, 0.1);
        sigmaSleep[i] ~ cauchy(0,2.5);
        sigmaActive[i] ~ cauchy(0,2.5);
        for(j in 2 : t ){
            if (Act[i,j]==0){
                target+=gamma_lpdf(nearZero|shape[i,j],rate[i,j]);
            }else{
                target+=gamma_lpdf(Act[i,j]|shape[i,j],rate[i,j]);
            }
            if (pEAP[i,j]<0.5){
                gammaMu[i,j] ~ normal(gammaMu[i,j-1],sigmaSleep[i]);
            }else{
                gammaMu[i,j] ~ normal(gammaMu[i,j-1],sigmaActive[i]);
        }
    }
    }
}

generated quantities{
    matrix <lower=0,upper=1> [N,t]predZeroOne;
    matrix <lower=0> [N,t]predAct;
    for(i in 1: N){
        for(j in 1 : t){
            predZeroOne[i,j]=bernoulli_rng(pEAP[i,j]);
            if(predZeroOne[i,j]==0){
                predAct[i,j]=0;
            }else{
                predAct[i,j]=gamma_rng(shape[i,j],rate[i,j]);
            }
        }
    }
}

"

datastan = list(N=N,t=t,Act=Act,pEAP=pEAP)

#fit
fit5 = stan(model_code = model_ZIG,
            data = datastan,
            chains = 4,
            iter = 20000,
            warmup = 4000,
            thin=32,
            seed = 1234)


#Apply Summarize Function
res6=mySummary(rstan::extract(fit5)$predAct)

#posterior prediction
ggplot(res6,aes(x=TimePoint,colour=ID,fill=ID))+
    geom_line(aes(y=EAP))+
    geom_line(aes(y=Activity))+
    geom_errorbar(aes(ymax=Upper95CI,ymin=Lower95CI),alpha=0.3)+
    facet_wrap(~ID,ncol = 1)+
    scale_x_continuous(expand = c(0,0))

7.3. Results

Figure is showed the posterior predictive distributions(errorbars indicate 95% Intervals).

Posterior predictive distributions illustrate the different activity variance between estimated sleep time and awake time.

Resut_ZIG_seed1234.png
Figure Zero-Inflated gamma dynamic model for human activity data.

8. Discussion

The normal state space model which is assumed gaussian system and gaussian observation could not describe human activity data, especially in successive zero-values. Non-gaussian dynamic models will performe better than normal model. This article proposed some possibility of Bayesian modeling for human activity using non-gaussian distribution. Berunoulli-Beta model and Zero-Inflated Gamma model described well the “sleep-awake probability” and the activity(and its state). However, our models yet including theoretical unknowns. Further studies are needed.

Written on July 19, 2017