Reports: Religious Change is a Simpson's Paradox

Joseph Bulbulia; Chris G. Sibley

Religious change at a societial level indicates “clocklike” decline

Figure ?? a. presents New Zealand Census data for society-level religious affiliation. The census data show a clock-like decline in religious affiliation at the level of the New Zealand population since 1966 of about 1% loss in religious affiliation per year, and more recently, closer to 2% loss. This suggests that people are becoming less religious over time.

Religious change at the individual level consists of opposing processes

However, society-level change is the average of two opposing processes that occur at the level of individuals:

religious affiliation: people might affiliate from a position of no-affiliation.
religious disaffiliation: people might disaffiliate from a position of affiliation.

These two opposing processes occur against a background of stability - most people remain constant in their religious affiliation or disaffiliation.

Here, we use a Hidden Markov Model to estimate rates of religious affiliation and disaffiliation across the New Zealand population from 2009/10 – 2020 in a sample of respondents who completed at least two NZAVS waves (\(N\) = 48272).

We find that, annually, the society-level rate of religious disaffiliation in New Zealand is .02%. That is the difference between the disaffiliation rate (Pr = .03) and the religious affiliation rate (Pr = .01). NZAVS data do not capture fluctuations from migration, hence, the NZAVS closely approximates the national census estimate for net religious change at the society level.

Show code

# Using markov models, show that change within people as they age differs from the society-wide trend
df <- readRDS(here::here("data", "df"))

length(unique(df$Id))

[1] 67690

Show code



#create dataframe
dat  <- df %>%
  dplyr::filter(YearMeasured == 1) %>% # filtering only people who have responded in a year
  dplyr::select(Wave, Age, Wave, Religious1 , Id, Relid, years) %>%
  dplyr::mutate(yearW  = as.numeric(Wave) - 1) %>%
  # filter(!is.na(Religious1))%>%
  group_by(Id) %>% filter(n() > 1) %>% # have responded to at least 2 x waves
  filter(n() != 0) %>%
  filter(n() != 1) %>%
  ungroup(Id) %>%
  group_by(Wave) %>%
  filter(!is.na(Religious1)) %>%
  ungroup() %>%
  group_by(Id) %>% filter(n() > 1) %>% # have responded to at least 2 x waves
  ungroup() %>%
  arrange(Id, yearW)
#nrow(df)

# unique ids
length(unique(dat$Id)) #46681

[1] 46681

Participants by waves

Show code

count_waves_participants <- function(x){
  out<-dplyr::count(tally(group_by(x, Id), sort = TRUE, name="number_waves"), number_waves)
  print(out)
}

count_waves_participants(dat)

# A tibble: 10 × 2
   number_waves     n
          <int> <int>
 1            2 23638
 2            3  3982
 3            4  5365
 4            5  2018
 5            6  2437
 6            7  3494
 7            8  2195
 8            9  1582
 9           10   680
10           11  1290

MSM Analysis

Show code

### MSM ANALYIS BELOW ###
# commenting out model to save time (which I've saved in a file, below)
# initial matrix for transitions
library("msm")

Q0 <- rbind(c(.9, .1),
            c(.1, .9))

# intensity matrix
Q0

# Initialise
crudeinR <-
  msm::crudeinits.msm(Religious1 ~ yearW, Id, # Id = identifier
                      data = dat,
                      qmatrix = Q0)

m1  <- msm(
  Religious1 ~ yearW,
  Id,
  data = dat,
  covariates =  ~ Age,
  qmatrix = crudeinR,
  ematrix = rbind(c(.1, .1), c(.1, .1)),
  est.initprobs = TRUE,
  exacttimes = TRUE,
  # Not running without this assumption.
  gen.inits = TRUE
)


#saveRDS(m1, here::here( "models","mnm-simpson-paradox-religious-change.rds" ))
saveRDS(m1, here::here( "_posts","religious_simpsons_paradox", "models","mnm-simpson-paradox-religious-change.rds" ))

Transition probabilities (averaged)

Show code

# read model, (saved to save time)
m1 <-readRDS(here::here( "_posts","religious_simpsons_paradox", "models","mnm-simpson-paradox-religious-change.rds" ))
m1


Call:
msm(formula = Religious1 ~ yearW, subject = Id, data = dat, qmatrix = crudeinR,     gen.inits = TRUE, ematrix = rbind(c(0.1, 0.1), c(0.1, 0.1)),     covariates = ~Age, est.initprobs = TRUE, exacttimes = TRUE)

Maximum likelihood estimates
Baselines are with covariates set to their means

Transition intensities with hazard ratios for each covariate
                  Baseline                       
State 1 - State 1 -0.009038 (-0.009986,-0.008181)
State 1 - State 2  0.009038 ( 0.008181, 0.009986)
State 2 - State 1  0.018906 ( 0.017204, 0.020776)
State 2 - State 2 -0.018906 (-0.020776,-0.017204)
                  Age                   
State 1 - State 1                       
State 1 - State 2 1.0417 (1.0351,1.0484)
State 2 - State 1 0.9876 (0.9816,0.9937)
State 2 - State 2                       

Misclassification probabilities
                      Baseline                 
Obs State 1 | State 1 0.95800 (0.95625,0.95968)
Obs State 2 | State 1 0.04200 (0.04032,0.04375)
Obs State 1 | State 2 0.06214 (0.05939,0.06501)
Obs State 2 | State 2 0.93786 (0.93499,0.94061)

Initial state occupancy probabilities
         Estimate       LCL       UCL
State 1 0.6328347 0.6278637 0.6377438
State 2 0.3671653 0.3622562 0.3721363

-2 * log-likelihood:  140700.7

Show code

#probability matrix
dt<-  print( round( msm::pmatrix.msm( m1 ) , 2 ) )

        State 1 State 2
State 1    0.99    0.01
State 2    0.02    0.98

Show code


# label cols, rows:

rownames(dt) = c("Not-religious time 0", "Religious time 0")
colnames(dt) = c("Not-religious t+1","Religious t+1")

#table
dt%>%
 kbl(caption = "Hidden Markov Model: inferred population average religious state change, years 2009-2020") %>%
  kable_material(c("hover"))

Table 1: Hidden Markov Model: inferred population average religious state change, years 2009-2020
	Not-religious t+1	Religious t+1
Not-religious time 0	0.99	0.01
Religious time 0	0.02	0.98

LaTeX table

Show code

dt %>%
  kbl(booktabs = T,
      format = "latex",
      caption = "Hidden Markov Model: inferred population average religious state change, years 2009-2020")

We can compare this with the raw (un-modelled) state transition, which overstates overall change.

Show code

# example change in wave 2018/ wave 2019
dat2 <- dat %>%
  filter(yearW == 09 | yearW == 10)

raw <-msm::statetable.msm(Religious1, Id, data = dat2)
# calculates proportions by rows
out<- round( prop.table(raw, 1), 2)

rownames(out) = c("Not-religious time 0", "Religious time 0")
colnames(out) = c("Not-religious t+1","Religious t+1")

out %>%
 kbl(caption = "Raw sample religious state change: years 2009-2020") %>%
  kable_material(c("hover"))

Table 2: Raw sample religious state change: years 2009-2020
	Not-religious t+1	Religious t+1
Not-religious time 0	0.94	0.06
Religious time 0	0.16	0.84

LaTeX table

Show code

out %>%
  kbl(booktabs = T,
      format = "latex",
      caption ="Raw sample religious state change: years 2009-2020")



<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'># latex</span>
<span class='fu'><a href='https://rdrr.io/pkg/xtable/man/xtable.html'>xtable</a></span><span class='op'>(</span><span class='va'>out</span><span class='op'>)</span>
</code></pre></div>

</details>

</div>




<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'># Graph</span>
<span class='va'>oz</span><span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/matrix.html'>as.matrix</a></span><span class='op'>(</span> <span class='fu'>msm</span><span class='fu'>::</span><span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span> <span class='va'>m1</span> <span class='op'>)</span> <span class='op'>)</span>
<span class='va'>oz</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span><span class='va'>oz</span>, digits <span class='op'>=</span> <span class='fl'>2</span><span class='op'>)</span>


<span class='va'>stateNames</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"Secular (.63)"</span>,<span class='st'>"Religious (.37)"</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/base/row.names.html'>row.names</a></span><span class='op'>(</span><span class='va'>oz</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>; <span class='fu'><a href='https://rdrr.io/r/base/colnames.html'>colnames</a></span><span class='op'>(</span><span class='va'>oz</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>
<span class='va'>oz</span>
</code></pre></div>

</details>

            Secular (.63) Religious (.37)

Secular (.63) 0.99 0.01 Religious (.37) 0.02 0.98


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'>## need to transpose or else arrows go in the wrong direction</span>
<span class='va'>toz</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/t.html'>t</a></span><span class='op'>(</span><span class='va'>oz</span><span class='op'>)</span> 

<span class='co'>#library(pracma)</span>
<span class='fu'><a href='https://rdrr.io/pkg/diagram/man/plotmat.html'>plotmat</a></span><span class='op'>(</span><span class='va'>toz</span>,
       <span class='co'># pos = c(1,2), </span>
       <span class='co'># lwd = 1, </span>
       <span class='co'># box.lwd = 2, </span>
        cex.txt <span class='op'>=</span> <span class='fl'>0.9</span>, 
        box.size <span class='op'>=</span> <span class='fl'>0.09</span>, 
       <span class='co'># box.type = "circle", </span>
        box.prop <span class='op'>=</span> <span class='fl'>0.5</span>,
        arr.pos <span class='op'>=</span>  <span class='fl'>0.4</span>, 
        shadow.size <span class='op'>=</span> <span class='fl'>0.001</span>,
        box.col <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"lightgreen"</span>,<span class='st'>"lightblue"</span>,<span class='st'>"orange"</span><span class='op'>)</span>,
       <span class='co'># arr.length=.1,</span>
       <span class='co'># arr.width=.1,</span>
        self.cex <span class='op'>=</span> <span class='fl'>.4</span>,
        self.shifty <span class='op'>=</span> <span class='op'>-</span><span class='fl'>.01</span>,
        self.shiftx <span class='op'>=</span> <span class='fl'>.12</span>,
        main <span class='op'>=</span> <span class='st'>"Hidden Markov Model reveals religious switching
       New Zealand: years 2009-2020"</span><span class='op'>)</span>
</code></pre></div>

</details><div class="figure">
<img src="simpsons-paradox_files/figure-html5/hmmchange-1.png" alt="Hidden Markov Model reveals religious change within people differs from the society wide-aggregate trend. As people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from greater disaffiliation among young people than affiliation." width="1536" />
<p class="caption">(\#fig:hmmchange)Hidden Markov Model reveals religious change within people differs from the society wide-aggregate trend. As people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from greater disaffiliation among young people than affiliation.</p>
</div>

</div>



## Graph 


<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='co'># note to self, this code can be improved</span>
<span class='co'>### CREATE PLOT  ###  CONVERSION AGE PLOT</span>
<span class='co'>## 20 year olds</span>

<span class='co'>#m1 &lt;-readRDS(here::here( "models","mnm-simpson-paradox-religious-change" ))</span>

<span class='va'>p20</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>20</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p20.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>p20.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p20.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.20</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>20</span>, <span class='va'>p20.Convert</span>, <span class='va'>p20.Convert.L</span>, <span class='va'>p20.Convert.U</span><span class='op'>)</span>
<span class='co'>## 30 year olds</span>
<span class='va'>p30</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>30</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p30.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p30.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p30.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'>#combine</span>
<span class='va'>con.30</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>30</span>, <span class='va'>p30.Convert</span>, <span class='va'>p30.Convert.L</span>, <span class='va'>p30.Convert.U</span><span class='op'>)</span>
<span class='co'>## 40 year old</span>
<span class='va'>p40</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>40</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p40.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p40.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p40.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.40</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>40</span>, <span class='va'>p40.Convert</span>, <span class='va'>p40.Convert.L</span>, <span class='va'>p40.Convert.U</span><span class='op'>)</span>
<span class='co'>## 50 year old</span>
<span class='va'>p50</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>50</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p50.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p50.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p50.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.50</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>50</span>, <span class='va'>p50.Convert</span>, <span class='va'>p50.Convert.L</span>, <span class='va'>p50.Convert.U</span><span class='op'>)</span>
<span class='co'>## 60 year olds</span>
<span class='va'>p60</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>60</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p60.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p60.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p60.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.60</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>60</span>, <span class='va'>p60.Convert</span>, <span class='va'>p60.Convert.L</span>, <span class='va'>p60.Convert.U</span><span class='op'>)</span>
<span class='co'>## 70 year olds</span>
<span class='va'>p70</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>70</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p70.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p70.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p70.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.70</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>70</span>, <span class='va'>p70.Convert</span>, <span class='va'>p70.Convert.L</span>, <span class='va'>p70.Convert.U</span><span class='op'>)</span>
<span class='co'>### ## 80 year olds</span>
<span class='va'>p80</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>80</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p80.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># set lower bound</span>
<span class='va'>p80.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># get upper bound</span>
<span class='va'>p80.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.80</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>80</span>, <span class='va'>p80.Convert</span>, <span class='va'>p80.Convert.L</span>, <span class='va'>p80.Convert.U</span><span class='op'>)</span>
<span class='co'># ### ## 90 year olds</span>
<span class='va'>p90</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>90</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># # Get Transition to Relid Estimate</span>
<span class='va'>p90.Convert</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='co'># # set lower bound</span>
<span class='va'>p90.Convert.L</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>p90.Convert.U</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>1</span>, <span class='fl'>2</span><span class='op'>]</span>
<span class='va'>con.90</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>90</span>, <span class='va'>p90.Convert</span>, <span class='va'>p90.Convert.L</span>, <span class='va'>p90.Convert.U</span><span class='op'>)</span>
<span class='co'># ### create data frame for conversion plot</span>
<span class='va'>conplot</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/r/base/data.frame.html'>data.frame</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/cbind.html'>rbind</a></span><span class='op'>(</span><span class='va'>con.20</span>, <span class='va'>con.30</span>, <span class='va'>con.40</span>, <span class='va'>con.50</span>, <span class='va'>con.60</span>, <span class='va'>con.70</span>, <span class='va'>con.80</span>, <span class='va'>con.90</span><span class='op'>)</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/base/colnames.html'>colnames</a></span><span class='op'>(</span><span class='va'>conplot</span><span class='op'>)</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"age"</span>, <span class='st'>"conv_probability"</span>, <span class='st'>"conv_lower"</span>, <span class='st'>"conv_upper"</span><span class='op'>)</span>
<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://ggplot2.tidyverse.org'>ggplot2</a></span><span class='op'>)</span>
<span class='co'>#library(wesanderson)</span>
<span class='co'>#wes_palette("Zissou1", type = c("discrete"))</span>
<span class='co'>#library("ggsci")</span>

<span class='va'>plot.con</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span><span class='op'>(</span><span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span><span class='op'>(</span>
    x <span class='op'>=</span> <span class='va'>age</span>,
    y <span class='op'>=</span> <span class='va'>conv_probability</span>,
    ymin <span class='op'>=</span> <span class='va'>conv_lower</span>,
    ymax <span class='op'>=</span> <span class='va'>conv_upper</span>
  <span class='op'>)</span>,
  data <span class='op'>=</span> <span class='va'>conplot</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_point.html'>geom_point</a></span><span class='op'>(</span>stroke <span class='op'>=</span> <span class='fl'>2</span>, colour <span class='op'>=</span> <span class='st'>"#d1495b"</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_linerange.html'>geom_linerange</a></span><span class='op'>(</span>alpha <span class='op'>=</span>   <span class='fl'>1</span>, colour <span class='op'>=</span> <span class='st'>"#d1495b"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_smooth.html'>geom_smooth</a></span><span class='op'>(</span>colour <span class='op'>=</span> <span class='st'>"#d1495b"</span>,
              method <span class='op'>=</span> <span class='st'>"loess"</span>,
              se <span class='op'>=</span> <span class='cn'>FALSE</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_continuous.html'>scale_y_continuous</a></span><span class='op'>(</span>limits <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>0</span>, <span class='fl'>.1</span><span class='op'>)</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>x <span class='op'>=</span> <span class='st'>"Age"</span>, y <span class='op'>=</span> <span class='st'>"Religious Affiliation Transition Probability w/ 95% CI"</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggtheme.html'>theme_bw</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>"Religious (Dis)affiliation Dynamics in New Zealand: years 2009-2020"</span><span class='op'>)</span>

<span class='co'>#plot.con</span>
<span class='co'>#m1</span>

<span class='co'>##### CREAT PLOT FOR DECONVERSION</span>
<span class='co'>## 20 year olds</span>
<span class='va'>p20</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>20</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='co'># Get Transition to Relid Estimate</span>
<span class='va'>p20.decon</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p20.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p20.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p20</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.20</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>20</span>, <span class='va'>p20.decon</span>, <span class='va'>p20.decon.L</span>, <span class='va'>p20.decon.U</span><span class='op'>)</span>
<span class='va'>p30</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>30</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p30.decon</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p30.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p30.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p30</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.30</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>30</span>, <span class='va'>p30.decon</span>, <span class='va'>p30.decon.L</span>, <span class='va'>p30.decon.U</span><span class='op'>)</span>
<span class='va'>p40</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>40</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p40.decon</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p40.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p40.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p40</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.40</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>40</span>, <span class='va'>p40.decon</span>, <span class='va'>p40.decon.L</span>, <span class='va'>p40.decon.U</span><span class='op'>)</span>
<span class='va'>p50</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>50</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p50.decon</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p50.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p50.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p50</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.50</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>50</span>, <span class='va'>p50.decon</span>, <span class='va'>p50.decon.L</span>, <span class='va'>p50.decon.U</span><span class='op'>)</span>
<span class='va'>p60</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>60</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p60.decon</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p60.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p60.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p60</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.60</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>60</span>, <span class='va'>p60.decon</span>, <span class='va'>p60.decon.L</span>, <span class='va'>p60.decon.U</span><span class='op'>)</span>
<span class='va'>p70</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>70</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p70.decon</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p70.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p70.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p70</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.70</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>70</span>, <span class='va'>p70.decon</span>, <span class='va'>p70.decon.L</span>, <span class='va'>p70.decon.U</span><span class='op'>)</span>
<span class='va'>p80</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>80</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p80.decon</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p80.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p80.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p80</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.80</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>80</span>, <span class='va'>p80.decon</span>, <span class='va'>p80.decon.L</span>, <span class='va'>p80.decon.U</span><span class='op'>)</span>
<span class='va'>p90</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>90</span><span class='op'>)</span>, ci <span class='op'>=</span> <span class='st'>"normal"</span><span class='op'>)</span>
<span class='va'>p90.decon</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>estimates</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p90.decon.L</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>L</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>p90.decon.U</span> <span class='op'>&lt;-</span> <span class='va'>p90</span><span class='op'>$</span><span class='va'>U</span><span class='op'>[</span><span class='fl'>2</span>, <span class='fl'>1</span><span class='op'>]</span>
<span class='va'>decon.90</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>90</span>, <span class='va'>p90.decon</span>, <span class='va'>p90.decon.L</span>, <span class='va'>p90.decon.U</span><span class='op'>)</span>

<span class='co'>### create data frame for conversion plot</span>
<span class='va'>deconplot</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/r/base/data.frame.html'>data.frame</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/cbind.html'>rbind</a></span><span class='op'>(</span>
    <span class='va'>decon.20</span>,
    <span class='va'>decon.30</span>,
    <span class='va'>decon.40</span>,
    <span class='va'>decon.50</span>,
    <span class='va'>decon.60</span>,
    <span class='va'>decon.70</span>,
    <span class='va'>decon.80</span>,
    <span class='va'>decon.90</span>
  <span class='op'>)</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/base/colnames.html'>colnames</a></span><span class='op'>(</span><span class='va'>deconplot</span><span class='op'>)</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"age"</span>, <span class='st'>"decon_probability"</span>, <span class='st'>"decon_lower"</span>, <span class='st'>"decon_upper"</span><span class='op'>)</span>
<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://ggplot2.tidyverse.org'>ggplot2</a></span><span class='op'>)</span>
<span class='va'>plot.decon</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span><span class='op'>(</span><span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span><span class='op'>(</span>
    x <span class='op'>=</span> <span class='va'>age</span>,
    y <span class='op'>=</span> <span class='va'>decon_probability</span>,
    ymin <span class='op'>=</span> <span class='va'>decon_lower</span>,
    ymax <span class='op'>=</span> <span class='va'>decon_upper</span>
  <span class='op'>)</span>,
  data <span class='op'>=</span> <span class='va'>deconplot</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_point.html'>geom_point</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span>,
             colour <span class='op'>=</span> <span class='st'>"#00798c"</span>,
             stroke <span class='op'>=</span> <span class='fl'>2</span><span class='op'>)</span> <span class='op'>+</span>   <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_linerange.html'>geom_linerange</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span>, colour <span class='op'>=</span> <span class='st'>"#00798c"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='co'>#scale_y_continuous(limits=c(0,.1))    +</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_smooth.html'>geom_smooth</a></span><span class='op'>(</span>colour <span class='op'>=</span> <span class='st'>"#00798c"</span>, se <span class='op'>=</span> <span class='cn'>T</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>x <span class='op'>=</span> <span class='st'>"Age"</span>, y <span class='op'>=</span> <span class='st'>"Religious Disaffiliation Transition Probability w/ 95% CI"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggtheme.html'>theme_bw</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>"Religious (Dis)affiliation Dynamics"</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/theme.html'>theme</a></span><span class='op'>(</span>legend.position <span class='op'>=</span> <span class='st'>"bottom"</span><span class='op'>)</span>

<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'>gridExtra</span><span class='op'>)</span>

<span class='co'>#grid.arrange(plot.con,plot.decon, ncol=2)</span>
<span class='co'>### SUPERIMPOSE</span>
<span class='va'>deconplot</span><span class='op'>$</span><span class='va'>Group</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/factor.html'>as.factor</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/rep.html'>rep</a></span><span class='op'>(</span><span class='st'>"Disaffiliation"</span>, <span class='fu'><a href='https://rdrr.io/r/base/nrow.html'>nrow</a></span><span class='op'>(</span><span class='va'>deconplot</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>conplot</span><span class='op'>$</span><span class='va'>Group</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/factor.html'>as.factor</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/rep.html'>rep</a></span><span class='op'>(</span><span class='st'>"Affiliation"</span>, <span class='fu'><a href='https://rdrr.io/r/base/nrow.html'>nrow</a></span><span class='op'>(</span><span class='va'>conplot</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>deconplot2</span> <span class='op'>&lt;-</span> <span class='va'>deconplot</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/mutate.html'>mutate</a></span><span class='op'>(</span>
    conv_probability <span class='op'>=</span> <span class='va'>decon_probability</span>,
    conv_lower <span class='op'>=</span> <span class='va'>decon_lower</span>,
    conv_upper <span class='op'>=</span> <span class='va'>decon_upper</span>
  <span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/select.html'>select</a></span><span class='op'>(</span><span class='op'>-</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='va'>decon_probability</span>, <span class='va'>decon_upper</span>, <span class='va'>decon_lower</span><span class='op'>)</span><span class='op'>)</span>

<span class='va'>comboplot</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/cbind.html'>rbind</a></span><span class='op'>(</span><span class='va'>conplot</span>, <span class='va'>deconplot2</span><span class='op'>)</span>
<span class='co'>#comboplot</span>
<span class='co'># get initial probabilities for title label</span>
<span class='co'>#pmatrix.msm(m1)</span>

<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://ggplot2.tidyverse.org'>ggplot2</a></span><span class='op'>)</span>

<span class='co'># plot.combo &lt;-</span>
<span class='co'>#   ggplot(</span>
<span class='co'>#     aes(</span>
<span class='co'>#       x = age,</span>
<span class='co'>#       y = conv_probability,</span>
<span class='co'>#       ymin = conv_lower,</span>
<span class='co'>#       ymax = conv_upper,</span>
<span class='co'>#       colour = Group</span>
<span class='co'>#     ),</span>
<span class='co'>#     data = comboplot</span>
<span class='co'>#   ) +</span>
<span class='co'>#   geom_point(alpha = 1, stroke = 2) +</span>
<span class='co'>#   geom_linerange(alpha = 1) +</span>
<span class='co'>#   geom_smooth(se = F) +</span>
<span class='co'>#   scale_y_continuous(limits = c(0, .96))  +</span>
<span class='co'>#   labs(x = "Age", y = "Annual Probability of Religious (Dis)affiliation w/ 95 CI: range: Pr = 0 to Pr = 1 ") + theme_bw() + labs(title = "Religious (Dis)affiliation Dynamics in New Zealand: years 2020 -- 2020") +</span>
<span class='co'>#   #scale_colour_manual(values=c( "#E69F00", "#56B4E9")) +</span>
<span class='co'>#   scale_colour_manual(values = c("#d1495b", "#00798c")) +</span>
<span class='co'>#   theme(legend.position = "bottom")</span>
<span class='co'># plot.combo</span>


<span class='va'>plotindividual</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span><span class='op'>(</span>
    <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span><span class='op'>(</span>
      x <span class='op'>=</span> <span class='va'>age</span>,
      y <span class='op'>=</span> <span class='va'>conv_probability</span>,
      ymin <span class='op'>=</span> <span class='va'>conv_lower</span>,
      ymax <span class='op'>=</span> <span class='va'>conv_upper</span>,
      colour <span class='op'>=</span> <span class='va'>Group</span>
    <span class='op'>)</span>,
    data <span class='op'>=</span> <span class='va'>comboplot</span>
  <span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_point.html'>geom_point</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span>, stroke <span class='op'>=</span> <span class='fl'>2</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_linerange.html'>geom_linerange</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_smooth.html'>geom_smooth</a></span><span class='op'>(</span>se <span class='op'>=</span> <span class='cn'>F</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_continuous.html'>scale_y_continuous</a></span><span class='op'>(</span>limits <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>0</span>, <span class='fl'>1</span><span class='op'>)</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>x <span class='op'>=</span> <span class='st'>"Age"</span>, y <span class='op'>=</span> <span class='st'>"Annual Probability of Religious (Dis)affiliation w/ 95 CI "</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggtheme.html'>theme_bw</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>"Religious (Dis)affiliation Dynamics in New Zealand: years 2009-2020"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='co'>#scale_colour_manual(values=c( "#E69F00", "#56B4E9")) +</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_manual.html'>scale_colour_manual</a></span><span class='op'>(</span>values <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"#d1495b"</span>, <span class='st'>"#00798c"</span><span class='op'>)</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/theme.html'>theme</a></span><span class='op'>(</span>legend.position <span class='op'>=</span> <span class='st'>"bottom"</span><span class='op'>)</span>

<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://patchwork.data-imaginist.com'>patchwork</a></span><span class='op'>)</span>
<span class='va'>pr</span> <span class='op'>+</span> <span class='va'>plotindividual</span>  <span class='op'>+</span> 
  <span class='fu'><a href='https://patchwork.data-imaginist.com/reference/plot_annotation.html'>plot_annotation</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>'Aggregated society-level change (a) differs from Individual-level religious change (b)'</span>, tag_levels <span class='op'>=</span> <span class='st'>"a"</span><span class='op'>)</span>
</code></pre></div>

</details><div class="figure">
<img src="simpsons-paradox_files/figure-html5/hmmchange2-1.png" alt="Hidden Markov Model reveals religious change within people differs from the society wide-aggregate trend. As people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation." width="1536" />
<p class="caption">(\#fig:hmmchange2)Hidden Markov Model reveals religious change within people differs from the society wide-aggregate trend. As people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation.</p>
</div><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">

<span class='co'>#pr / plot.combo + plot_annotation(title = "Religious Change over Time",</span>
<span class='co'>#         subtitle = "(i) Society-level aggregate change (ii) Individual-level process: affiliation vs. disaffiliation",</span>
<span class='co'>#tag_levels = "i")</span>
</code></pre></div>

</details>

</div>


However, as evident in Figure \@ref(fig:hmmchange),
individual-level NZAVS data reveal that as people age, the
religious affiliation rate gradually increases. Likewise, as
people age, the religious disaffiliation rate gradually
decreases. The age at which religious affiliation exceeds
disaffiliation is about 70 years old. The society-level
trend and individual-level trends run in the opposite
directions because there is higher relative religious
disaffiliation among young people. This example illustrates
what is called a "Simpson's Paradox." Wherever there is
potential for switching between states at an individual
level, there is on opportunity for a Simpson's Paradox. Most
facts of interest about the role of religion in a persons
life are subject to switching.

## Religious change is a Simpson's Paradox

<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='va'>plotI</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span><span class='op'>(</span>
    <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span><span class='op'>(</span>
      x <span class='op'>=</span> <span class='va'>age</span>,
      y <span class='op'>=</span> <span class='va'>conv_probability</span>,
      ymin <span class='op'>=</span> <span class='va'>conv_lower</span>,
      ymax <span class='op'>=</span> <span class='va'>conv_upper</span>,
      colour <span class='op'>=</span> <span class='va'>Group</span>
    <span class='op'>)</span>,
    data <span class='op'>=</span> <span class='va'>comboplot</span>
  <span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_point.html'>geom_point</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span>, stroke <span class='op'>=</span> <span class='fl'>2</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_linerange.html'>geom_linerange</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_smooth.html'>geom_smooth</a></span><span class='op'>(</span>se <span class='op'>=</span> <span class='cn'>F</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_continuous.html'>scale_y_continuous</a></span><span class='op'>(</span>limits <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>0</span>, <span class='fl'>.96</span><span class='op'>)</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>x <span class='op'>=</span> <span class='st'>"Age"</span>, y <span class='op'>=</span> <span class='st'>"Annual Probability of Religious (Dis)affiliation w/ 95 CI: range: Pr = 0 to Pr = 1 "</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggtheme.html'>theme_bw</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>"Religious (Dis)affiliation Dynamics in New Zealand: years 2009--2020"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='co'>#scale_colour_manual(values=c( "#E69F00", "#56B4E9")) +</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_manual.html'>scale_colour_manual</a></span><span class='op'>(</span>values <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"#d1495b"</span>, <span class='st'>"#00798c"</span><span class='op'>)</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/theme.html'>theme</a></span><span class='op'>(</span>legend.position <span class='op'>=</span> <span class='st'>"bottom"</span><span class='op'>)</span>
<span class='va'>plotI</span>
</code></pre></div>

</details><div class="figure">
<img src="simpsons-paradox_files/figure-html5/hmmchangeI-1.png" alt="Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation." width="1536" />
<p class="caption">(\#fig:hmmchangeI-1)Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation.</p>
</div><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='va'>plotIb</span> <span class='op'>&lt;-</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggplot.html'>ggplot</a></span><span class='op'>(</span>
    <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/aes.html'>aes</a></span><span class='op'>(</span>
      x <span class='op'>=</span> <span class='va'>age</span>,
      y <span class='op'>=</span> <span class='va'>conv_probability</span>,
      ymin <span class='op'>=</span> <span class='va'>conv_lower</span>,
      ymax <span class='op'>=</span> <span class='va'>conv_upper</span>,
      colour <span class='op'>=</span> <span class='va'>Group</span>
    <span class='op'>)</span>,
    data <span class='op'>=</span> <span class='va'>comboplot</span>
  <span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_point.html'>geom_point</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span>, stroke <span class='op'>=</span> <span class='fl'>2</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_linerange.html'>geom_linerange</a></span><span class='op'>(</span>alpha <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/geom_smooth.html'>geom_smooth</a></span><span class='op'>(</span>se <span class='op'>=</span> <span class='cn'>F</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_continuous.html'>scale_y_continuous</a></span><span class='op'>(</span>limits <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>0</span>, <span class='fl'>.1</span><span class='op'>)</span><span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>x <span class='op'>=</span> <span class='st'>"Age"</span>, y <span class='op'>=</span> <span class='st'>"Annual Probability of Religious (Dis)affiliation w/ 95 CI: range: Pr = 0 to Pr = .1 "</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/ggtheme.html'>theme_bw</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>+</span> <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/labs.html'>labs</a></span><span class='op'>(</span>title <span class='op'>=</span> <span class='st'>"Religious (Dis)affiliation Dynamics in New Zealand: years 2009--2020"</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='co'>#scale_colour_manual(values=c( "#E69F00", "#56B4E9")) +</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/scale_manual.html'>scale_colour_manual</a></span><span class='op'>(</span>values <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"#d1495b"</span>, <span class='st'>"#00798c"</span><span class='op'>)</span><span class='op'>)</span> <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/theme.html'>theme</a></span><span class='op'>(</span>legend.position <span class='op'>=</span> <span class='st'>"bottom"</span><span class='op'>)</span>

<span class='va'>plotIb</span>
</code></pre></div>

</details><div class="figure">
<img src="simpsons-paradox_files/figure-html5/hmmchangeI-2.png" alt="Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation." width="1536" />
<p class="caption">(\#fig:hmmchangeI-2)Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation.</p>
</div><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='va'>plotI</span> <span class='op'>/</span><span class='va'>plotIb</span>
</code></pre></div>

</details><div class="figure">
<img src="simpsons-paradox_files/figure-html5/hmmchangeI-3.png" alt="Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation." width="1536" />
<p class="caption">(\#fig:hmmchangeI-3)Close up of opposing processes within indiviuals (scale pr = 0 to pr = .1). Hidden Markov Model reveals that as people age, they are more likely to affiliate and less likely to disaffiliate. The population gap comes from much greater disaffiliation among young people than affiliation.</p>
</div><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'># # to save file in your directory at correct resolution, uncomment below</span>
<span class='co'>#  file &lt;- here::here("_posts","religious_simpsons_paradox", "figs", paste("aff_disaff", ".png", sep = ""))</span>
<span class='co'># #</span>
<span class='co'># ggsave(</span>
<span class='co'>#   file,</span>
<span class='co'>#   device = "png",</span>
<span class='co'>#   scale = 1,</span>
<span class='co'>#   width = 16,</span>
<span class='co'>#   height = 9,</span>
<span class='co'>#   units = "in",</span>
<span class='co'>#   dpi = 300,</span>
<span class='co'>#   limitsize = TRUE</span>
<span class='co'># )</span>
</code></pre></div>

</details>

</div>


Appreciating that religious change is a Simpson's paradox
helps to build an intuition for why panel data are so
powerful. Sociologists have documented persistent declines
in religious affiliation across western societies, a process
Max Weber called "the disenchantment of the world"
[@Weber1946-zp]. This disenchantment theory is a central
dogma in the sociology of religion [@Norris2004-gu]. The
theory is accurate for describing society-level dynamics.
Indeed, NZAVS data confirm the society-level
"disenchantment." However, this central dogma is poorly
suited for understanding how religion and spirituality take
shape within the lives of individuals.


## Appendix

## Descriptive Table



## Raw Trajectories

<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'>LMest</span><span class='op'>)</span>
<span class='va'>relwidelong</span> <span class='op'>&lt;-</span>  <span class='va'>df</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/mutate.html'>mutate</a></span><span class='op'>(</span>Religious <span class='op'>=</span> <span class='va'>Religious1</span> <span class='op'>-</span> <span class='fl'>1</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/select.html'>select</a></span><span class='op'>(</span><span class='va'>Id</span>, <span class='va'>Religious</span>, <span class='va'>Wave</span><span class='op'>)</span>

<span class='co'>##THIS HANGS</span>
<span class='co'>#out &lt;- with(relwidelong, long2matrices(id = Id, X = Age, Y = Religious))</span>

<span class='va'>relwide</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://tidyr.tidyverse.org/reference/spread.html'>spread</a></span><span class='op'>(</span><span class='va'>relwidelong</span>, <span class='va'>Wave</span>, <span class='va'>Religious</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/utils/head.html'>head</a></span><span class='op'>(</span><span class='va'>relwide</span><span class='op'>)</span>
</code></pre></div>

</details>

A tibble: 6 × 12

Id 2009 2010 2011 2012 2013 2014 2015 1 1 0 NA NA NA NA NA NA 2 2 1 NA NA NA 0 0 0 3 3 0 NA NA NA NA NA NA 4 4 0 NA NA NA NA NA NA 5 5 1 0 NA NA NA NA NA 6 6 0 NA NA NA NA NA NA # … with 4 more variables: 2016 , 2017 , # 2018 , 2019


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='co'>#relwide&lt;-relwide[complete.cases(relwide), ]</span>


<span class='va'>Relidwide</span> <span class='op'>&lt;-</span>  <span class='va'>df</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span><span class='op'>(</span><span class='va'>YearMeasured</span> <span class='op'>==</span> <span class='fl'>1</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span><span class='op'>(</span><span class='op'>!</span><span class='fu'><a href='https://rdrr.io/r/base/NA.html'>is.na</a></span><span class='op'>(</span><span class='va'>Religious</span><span class='op'>)</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'>dplyr</span><span class='fu'>::</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/select.html'>select</a></span><span class='op'>(</span><span class='va'>Id</span>, <span class='va'>Religious</span>, <span class='va'>Wave</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'><a href='https://dplyr.tidyverse.org/reference/group_by.html'>group_by</a></span><span class='op'>(</span><span class='va'>Id</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span> <span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span><span class='op'>(</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/context.html'>n</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>&gt;</span> <span class='fl'>10</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'><a href='https://dplyr.tidyverse.org/reference/filter.html'>filter</a></span><span class='op'>(</span><span class='fu'><a href='https://dplyr.tidyverse.org/reference/context.html'>n</a></span><span class='op'>(</span><span class='op'>)</span> <span class='op'>!=</span> <span class='fl'>0</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'><a href='https://dplyr.tidyverse.org/reference/group_by.html'>ungroup</a></span><span class='op'>(</span><span class='va'>Id</span><span class='op'>)</span> <span class='op'><a href='https://rdrr.io/pkg/igraph/man/pipe.html'>%&gt;%</a></span>
  <span class='fu'><a href='https://tidyr.tidyverse.org/reference/spread.html'>spread</a></span><span class='op'>(</span><span class='va'>Wave</span>, <span class='va'>Religious</span><span class='op'>)</span>

<span class='va'>Relidwide</span> <span class='op'>&lt;-</span> <span class='va'>Relidwide</span><span class='op'>[</span><span class='fu'><a href='https://rdrr.io/r/stats/complete.cases.html'>complete.cases</a></span><span class='op'>(</span><span class='va'>Relidwide</span><span class='op'>)</span>,<span class='op'>]</span>
<span class='fu'><a href='https://rdrr.io/r/base/dim.html'>dim</a></span><span class='op'>(</span><span class='va'>Relidwide</span><span class='op'>)</span>
</code></pre></div>

</details>

[1] 1290 12


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='va'>x_var_list</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/names.html'>names</a></span><span class='op'>(</span><span class='va'>Relidwide</span><span class='op'>[</span>, <span class='fl'>2</span><span class='op'>:</span><span class='fl'>12</span><span class='op'>]</span><span class='op'>)</span>

<span class='va'>Relidwide</span>
</code></pre></div>

</details>

A tibble: 1,290 × 12

Id 2009 2010 2011 2012 2013 2014 2015 1 15 Not_Reli… Not_R… Not_R… Not_R… Not_R… Not_R… Not_R… 2 19 Not_Reli… Relig… Not_R… Relig… Relig… Relig… Relig… 3 20 Not_Reli… Not_R… Not_R… Not_R… Not_R… Not_R… Not_R… 4 28 Not_Reli… Not_R… Not_R… Not_R… Not_R… Not_R… Not_R… 5 34 Not_Reli… Relig… Not_R… Not_R… Relig… Relig… Relig… 6 38 Religious Relig… Relig… Relig… Relig… Relig… Relig… 7 39 Not_Reli… Not_R… Not_R… Not_R… Not_R… Not_R… Not_R… 8 40 Not_Reli… Not_R… Not_R… Not_R… Not_R… Not_R… Not_R… 9 41 Not_Reli… Not_R… Relig… Not_R… Not_R… Not_R… Not_R… 10 46 Religious Relig… Relig… Relig… Relig… Relig… Relig… # … with 1,280 more rows, and 4 more variables: # 2016 , 2017 , 2018 , 2019


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">
<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://milanwiedemann.github.io/lcsm/'>lcsm</a></span><span class='op'>)</span>


<span class='fu'><a href='https://milanwiedemann.github.io/lcsm/reference/plot_trajectories.html'>plot_trajectories</a></span><span class='op'>(</span>
  data <span class='op'>=</span> <span class='va'>Relidwide</span>,
  id_var <span class='op'>=</span> <span class='st'>"Id"</span>,
  var_list <span class='op'>=</span> <span class='va'>x_var_list</span>,
  line_colour <span class='op'>=</span> <span class='st'>"red"</span>,
  xlab <span class='op'>=</span> <span class='st'>"Time"</span>,
  ylab <span class='op'>=</span> <span class='st'>"Value"</span>,
  connect_missing <span class='op'>=</span> <span class='cn'>F</span>,
  scale_x_num <span class='op'>=</span> <span class='cn'>T</span>,
  scale_x_num_start <span class='op'>=</span> <span class='fl'>1</span>,
  random_sample_frac <span class='op'>=</span> <span class='fl'>0.05</span>,
  seed <span class='op'>=</span> <span class='fl'>1234</span>
 <span class='co'># title_n = T</span>
<span class='op'>)</span>  <span class='op'>+</span>
  <span class='fu'><a href='https://ggplot2.tidyverse.org/reference/facet_wrap.html'>facet_wrap</a></span><span class='op'>(</span> <span class='op'>~</span> <span class='va'>Id</span><span class='op'>)</span>
</code></pre></div>

</details><img src="simpsons-paradox_files/figure-html5/unnamed-chunk-7-1.png" width="1536" /><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">

<span class='co'># # to save file in your directory at correct resolution, uncomment below</span>
<span class='co'># file &lt;- here::here("_posts","religious_simpsons_paradox", "figs", paste("trajectories", ".png", sep = ""))</span>
<span class='co'># # #</span>
<span class='co'># ggsave(</span>
<span class='co'>#   file,</span>
<span class='co'>#   plot = last_plot(),</span>
<span class='co'>#   device = "png",</span>
<span class='co'>#   scale = 1,</span>
<span class='co'>#   width = 16,</span>
<span class='co'>#   height = 9,</span>
<span class='co'>#   units = "in",</span>
<span class='co'>#   dpi = 300,</span>
<span class='co'>#   limitsize = TRUE</span>
<span class='co'># )</span>
</code></pre></div>

</details>

</div>




### Graph of states 

<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'># Graph</span>
<span class='va'>oz_20</span><span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/matrix.html'>as.matrix</a></span><span class='op'>(</span> <span class='fu'>msm</span><span class='fu'>::</span><span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>20</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>oz_20</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span><span class='va'>oz_20</span>, digits <span class='op'>=</span> <span class='fl'>3</span><span class='op'>)</span>


<span class='va'>stateNames</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"Secular "</span>,<span class='st'>"Religious"</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/base/row.names.html'>row.names</a></span><span class='op'>(</span><span class='va'>oz_20</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>; <span class='fu'><a href='https://rdrr.io/r/base/colnames.html'>colnames</a></span><span class='op'>(</span><span class='va'>oz_20</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>
<span class='va'>oz_20</span>
</code></pre></div>

</details>

      Secular  Religious

Secular 0.997 0.003 Religious 0.027 0.973


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'>## need to transpose or else arrows go in the wrong direction</span>
<span class='va'>toz_20</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/t.html'>t</a></span><span class='op'>(</span><span class='va'>oz_20</span><span class='op'>)</span> 

<span class='co'>#library(pracma)</span>
<span class='fu'><a href='https://rdrr.io/pkg/diagram/man/plotmat.html'>plotmat</a></span><span class='op'>(</span><span class='va'>toz_20</span>,
       <span class='co'># pos = c(1,2), </span>
       <span class='co'># lwd = 1, </span>
       <span class='co'># box.lwd = 2, </span>
        cex.txt <span class='op'>=</span> <span class='fl'>0.9</span>, 
        box.size <span class='op'>=</span> <span class='fl'>0.09</span>, 
       <span class='co'># box.type = "circle", </span>
        box.prop <span class='op'>=</span> <span class='fl'>0.5</span>,
        arr.pos <span class='op'>=</span>  <span class='fl'>0.4</span>, 
        shadow.size <span class='op'>=</span> <span class='fl'>0.001</span>,
        box.col <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"lightgreen"</span>,<span class='st'>"lightblue"</span>,<span class='st'>"orange"</span><span class='op'>)</span>,
       <span class='co'># arr.length=.1,</span>
       <span class='co'># arr.width=.1,</span>
        self.cex <span class='op'>=</span> <span class='fl'>.4</span>,
        self.shifty <span class='op'>=</span> <span class='op'>-</span><span class='fl'>.01</span>,
        self.shiftx <span class='op'>=</span> <span class='fl'>.12</span>,
        main <span class='op'>=</span> <span class='st'>"Hidden Markov Model reveals religious switching
       New Zealand: years 2009-2020, Age = 20"</span><span class='op'>)</span>
</code></pre></div>

</details><img src="simpsons-paradox_files/figure-html5/unnamed-chunk-8-1.png" width="1536" /><details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">


<span class='co'># Graph</span>
<span class='va'>oz_70</span><span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/matrix.html'>as.matrix</a></span><span class='op'>(</span> <span class='fu'>msm</span><span class='fu'>::</span><span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m1</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>70</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>oz_70</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span><span class='va'>oz_70</span>, digits <span class='op'>=</span> <span class='fl'>3</span><span class='op'>)</span>


<span class='va'>stateNames</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"Secular "</span>,<span class='st'>"Religious"</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/r/base/row.names.html'>row.names</a></span><span class='op'>(</span><span class='va'>oz_70</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>; <span class='fu'><a href='https://rdrr.io/r/base/colnames.html'>colnames</a></span><span class='op'>(</span><span class='va'>oz_70</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='va'>stateNames</span>
<span class='va'>oz_70</span>
</code></pre></div>

</details>

      Secular  Religious

Secular 0.980 0.020 Religious 0.015 0.985


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='co'>## need to transpose or else arrows go in the wrong direction</span>
<span class='va'>toz_70</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/t.html'>t</a></span><span class='op'>(</span><span class='va'>oz_70</span><span class='op'>)</span> 

<span class='co'>#library(pracma)</span>
<span class='fu'><a href='https://rdrr.io/r/grDevices/dev.html'>dev.off</a></span><span class='op'>(</span><span class='op'>)</span>
</code></pre></div>

</details>

null device 1


<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='fu'><a href='https://rdrr.io/r/graphics/par.html'>par</a></span><span class='op'>(</span>mfrow<span class='op'>=</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>1</span>,<span class='fl'>2</span><span class='op'>)</span><span class='op'>)</span>

<span class='fu'><a href='https://rdrr.io/pkg/diagram/man/plotmat.html'>plotmat</a></span><span class='op'>(</span><span class='va'>toz_20</span>,
       <span class='co'># pos = c(1,2), </span>
       <span class='co'># lwd = 1, </span>
       <span class='co'># box.lwd = 2, </span>
        cex.txt <span class='op'>=</span> <span class='fl'>0.9</span>, 
        box.size <span class='op'>=</span> <span class='fl'>0.09</span>, 
       <span class='co'># box.type = "circle", </span>
        box.prop <span class='op'>=</span> <span class='fl'>0.5</span>,
        arr.pos <span class='op'>=</span>  <span class='fl'>0.4</span>, 
        shadow.size <span class='op'>=</span> <span class='fl'>0.001</span>,
        box.col <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"lightgreen"</span>,<span class='st'>"lightblue"</span>,<span class='st'>"orange"</span><span class='op'>)</span>,
       <span class='co'># arr.length=.1,</span>
       <span class='co'># arr.width=.1,</span>
        self.cex <span class='op'>=</span> <span class='fl'>.4</span>,
        self.shifty <span class='op'>=</span> <span class='op'>-</span><span class='fl'>.01</span>,
        self.shiftx <span class='op'>=</span> <span class='fl'>.12</span>,
        main <span class='op'>=</span> <span class='st'>"Hidden Markov Model reveals religious switching
       New Zealand: years 2009-2020, Age = 20"</span><span class='op'>)</span>


<span class='fu'><a href='https://rdrr.io/pkg/diagram/man/plotmat.html'>plotmat</a></span><span class='op'>(</span><span class='va'>toz_70</span>,
       <span class='co'># pos = c(1,2), </span>
       <span class='co'># lwd = 1, </span>
       <span class='co'># box.lwd = 2, </span>
        cex.txt <span class='op'>=</span> <span class='fl'>0.9</span>, 
        box.size <span class='op'>=</span> <span class='fl'>0.09</span>, 
       <span class='co'># box.type = "circle", </span>
        box.prop <span class='op'>=</span> <span class='fl'>0.5</span>,
        arr.pos <span class='op'>=</span>  <span class='fl'>0.4</span>, 
        shadow.size <span class='op'>=</span> <span class='fl'>0.001</span>,
        box.col <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"lightgreen"</span>,<span class='st'>"lightblue"</span>,<span class='st'>"orange"</span><span class='op'>)</span>,
       <span class='co'># arr.length=.1,</span>
       <span class='co'># arr.width=.1,</span>
        self.cex <span class='op'>=</span> <span class='fl'>.4</span>,
        self.shifty <span class='op'>=</span> <span class='op'>-</span><span class='fl'>.01</span>,
        self.shiftx <span class='op'>=</span> <span class='fl'>.12</span>,
        main <span class='op'>=</span> <span class='st'>"Hidden Markov Model reveals religious switching
       New Zealand: years 2009-2020, Age = 70"</span><span class='op'>)</span>

<span class='co'>#age_compare 16 x 9</span>

<span class='co'># </span>
<span class='co'># </span>
<span class='co'># file &lt;- here::here("_posts","religious_simpsons_paradox", "figs", paste("age20_age70", ".png", sep = ""))</span>
<span class='co'># ##</span>
<span class='co'># </span>
<span class='co'># ggsave(</span>
<span class='co'>#   file,</span>
<span class='co'>#   plot = last_plot(),</span>
<span class='co'>#   device = "png",</span>
<span class='co'>#   scale = 1,</span>
<span class='co'>#   width = 16,</span>
<span class='co'>#   height = 9,</span>
<span class='co'>#   units = "in",</span>
<span class='co'>#   dpi = 300,</span>
<span class='co'>#   limitsize = TRUE</span>
<span class='co'># )</span>
</code></pre></div>

</details>

</div>






### Three-state model 


<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='va'>q3</span> <span class='op'>&lt;-</span><span class='fu'><a href='https://rdrr.io/r/base/cbind.html'>rbind</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>,<span class='fl'>.1</span>,<span class='fl'>.1</span><span class='op'>)</span>,
           <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>,<span class='fl'>.1</span>,<span class='fl'>.1</span><span class='op'>)</span>,
           <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>,<span class='fl'>.1</span>,<span class='fl'>.1</span><span class='op'>)</span><span class='op'>)</span> <span class='co'>#  Allow permitted transitions</span>


<span class='co'>#three state emissions model</span>
<span class='va'>m3</span>  <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/msm/man/msm.html'>msm</a></span><span class='op'>(</span>
  <span class='va'>Religious1</span> <span class='op'>~</span> <span class='va'>yearW</span>,
  <span class='va'>Id</span>,
  data <span class='op'>=</span> <span class='va'>dat</span>,
  covariates <span class='op'>=</span> <span class='op'>~</span> <span class='va'>Age</span>,
  est.initprobs <span class='op'>=</span> <span class='cn'>TRUE</span>,
  exacttimes <span class='op'>=</span> <span class='cn'>TRUE</span>,<span class='co'># Not running without this assumption.</span>
  gen.inits <span class='op'>=</span> <span class='cn'>TRUE</span>,
  qmatrix <span class='op'>=</span> <span class='va'>q3</span>,
  ematrix <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/cbind.html'>rbind</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>, <span class='fl'>.1</span>, <span class='fl'>.1</span><span class='op'>)</span>,
                  <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>, <span class='fl'>.1</span>, <span class='fl'>.1</span><span class='op'>)</span>,
                  <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>.1</span>, <span class='fl'>.1</span>, <span class='fl'>.1</span><span class='op'>)</span>
<span class='op'>)</span><span class='op'>)</span>

<span class='co'># saved to save time</span>
<span class='co'>#saveRDS(m3, here::here("models","simpsons-paradox-religion-three-state"))</span>
</code></pre></div>

</details>

</div>



<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='va'>m3</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/readRDS.html'>readRDS</a></span><span class='op'>(</span><span class='fu'>here</span><span class='fu'>::</span><span class='fu'><a href='https://here.r-lib.org//reference/here.html'>here</a></span><span class='op'>(</span><span class='st'>"models"</span>, <span class='st'>"simpsons-paradox-religion-three-state"</span><span class='op'>)</span><span class='op'>)</span>

<span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span> <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m3</span><span class='op'>)</span>, <span class='fl'>2</span><span class='op'>)</span> <span class='co'># works</span>

<span class='va'>three_state_mean</span> <span class='op'>&lt;-</span><span class='fu'><a href='https://rdrr.io/r/base/data.frame.html'>data.frame</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span> <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m3</span>, covariates <span class='op'>=</span> <span class='st'>"mean"</span><span class='op'>)</span>, <span class='fl'>2</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>three_state_mean</span>
</code></pre></div>

</details>

</div>



<div class="layout-chunk" data-layout="l-body-outset">
<details>
<summary>Show code</summary>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class='va'>three_state_30</span> <span class='op'>&lt;-</span><span class='fu'><a href='https://rdrr.io/r/base/data.frame.html'>data.frame</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span> <span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m3</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>30</span><span class='op'>)</span><span class='op'>)</span>, <span class='fl'>2</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>
<span class='va'>three_state_50</span> <span class='op'>&lt;-</span><span class='fu'><a href='https://rdrr.io/r/base/data.frame.html'>data.frame</a></span><span class='op'>(</span> <span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span> <span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/pkg/msm/man/pmatrix.msm.html'>pmatrix.msm</a></span><span class='op'>(</span><span class='va'>m3</span>, covariates <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>Age <span class='op'>=</span> <span class='fl'>50</span><span class='op'>)</span><span class='op'>)</span>, <span class='fl'>2</span><span class='op'>)</span><span class='op'>)</span><span class='op'>)</span>

<span class='va'>three_state_30</span>
<span class='va'>three_state_50</span>
<span class='fu'><a href='https://rdrr.io/pkg/xtable/man/xtable.html'>xtable</a></span><span class='op'>(</span><span class='va'>three_state_30</span><span class='op'>)</span>
<span class='fu'><a href='https://rdrr.io/pkg/xtable/man/xtable.html'>xtable</a></span><span class='op'>(</span><span class='va'>three_state_50</span><span class='op'>)</span>

<span class='co'>#round( pmatrix.msm(m3, covariates = "mean" ), 2)</span>
</code></pre></div>

</details>

</div>


```{.r .distill-force-highlighting-css}

Religious Change is a Simpson’s Paradox

Religious change at a societial level indicates “clocklike” decline

Religious change at the individual level consists of opposing processes

A tibble: 6 × 12

A tibble: 1,290 × 12

References

Corrections

Reuse

Citation