===== Daily religious study using multilevel mediation models =====
==== Overall purpose ====


Johnny – you can proceed in whatever fashion you want, but, one thing I want you to look at is the following mediation models

Relmin->theisf->Outcome

Relq->Theisq->Outcome

Relc->Theisc->Outcome

**The outcomes of most interest are Empx, Empc, Gratx, Gratc, Forx, and Forc. But, at some point I would also like to look at Mix, Mic, Intx, Intc, Purpx, Purpc.
===== Models =====
==== Without moderator ====

    - M1: Relmin --> theisf --> Empx
    - M2: Relmin --> theisf --> Empc
    - M3: Relmin --> theisf --> Gratx
    - M4: Relmin --> theisf --> Gratc
    - M5: Relmin --> theisf --> Forx
    - M6: Relmin --> theisf --> Forc
    - M7: Relq --> Theisq --> Empx
    - M8: Relq --> Theisq --> Empc
    - M9: Relq --> Theisq --> Gratx
    - M10: Relq --> Theisq --> Gratc
    - M11: Relq --> Theisq --> Forx
    - M12: Relq --> Theisq --> Forc
    - M13: Relc --> Theisc --> Empx
    - M14: Relc --> Theisc --> Empc
    - M15: Relc --> Theisc --> Gratx
    - M16: Relc --> Theisc --> Gratc
    - M17: Relc --> Theisc --> Forx
    - M18: Relc --> Theisc --> Forc

==== Data process ====
  - First select the data from Relmin, theisf, empx for a mediation model without covariate for M1
<code>
## Mediation analysis data preparation

setwd('/Users/zzhang4/zzy/research/Daily religion study/data')

dset<-read.table('DMA-data-2010-01-25.txt', header=T)
N<-dim(dset)[1]
case<-1:N

## for Relmin, theisf, empx

dset1<-dset[, c(1, 24, 12, 15, 18)]
dset1<-cbind(case, dset1)

## check for missing data
m<-is.na(dset1)

m.sum<-apply(m, 1, sum)

## take out individual with at least two data points
dset2<-dset1[m.sum<2, ]

## count the number of individual with certain number of data
count<-table(dset2$id)

dset2$freq.id<-rep(count, count)

## only keep individual with at least 21 data points
dset3<-dset2[dset2$freq.id>20, ]

## total number of individuals in the sample
freqid<-table(dset3$id)
length(table(dset3$id))->J
NN<-dim(dset3)[1]
id<-rep(1:J, freqid)
## Save data for mediation model
R<-diag(5)
## load the function rtobugs.r
source("http://nd.psychstat.org/_media/classdata/rtobugs.r")

tobugs<-list(N=NN, J=J, R=R, x=round(log(dset3$relmins),5), m=round(log(dset3$theisfs+1),5), y=dset3$empxm, ind=id)

rtobugs(tobugs, 'm1.txt')

</code>
==== Analysis using WinBUGS ====

[[A working example codes]]