Version 1: Strength of memory encoding and rank of items affect delayed inference performance

Dataset 1 (Within-subject Transitive Inference with 27h and 3h delayed testing) by Tamas n = 70

best fitting model is the one that includes interaction between Hierarchy ( time of encoding ) and Immediate premise pair performance ( strength of memory encoding ). Interaction is significant
main effect of distance, but not interacting with time of encoding overall
- there is main effect of time of encoding on the distant inference pair when looked at separately
there is a 3way interaction between strength of memory encoding, time of encoding and joint rank of items hinting at a certain kind of “spatial” time-dependent effect on inference performance
- the time of encoding and joint rank of items two-way interaction is also sig.

exp1.facescene <- readRDS("bytrial.exp-wstionline_ver-facescene.rds")
exp1.scenes <- readRDS("bytrial.exp-wstionline_ver-scenes.rds")
exp1.objects1 <- readRDS("bytrial.exp-wstionline_ver-objects.rds")
exp1.objects2 <- readRDS("bytrial.exp-wstionline_ver-objects-shp.rds")

# note to self - fix feedback column in exp1.objects
exp1.combined <- bind_rows(exp1.facescene,exp1.scenes,exp1.objects1,exp1.objects2) %>% mutate(participant_u = paste(participant,expVersion, sep = "_"), cRank= paste0(pmin(Rank1,Rank2),pmax(Rank1,Rank2)), pairType= ifelse(cRank %in% c("12","23","34","45","56"),"premise", ifelse(cRank %in% c("24","25","35"), "inference","anchor"))) %>% filter(!participant_u %in% c("ADPOMZ_objects","AKOXYS_objects_shp") )

write.csv(exp1.combined,"ellrep_ws.csv")

glimpse(exp1.combined)

## Rows: 26,040
## Columns: 18
## Groups: participant, hierarchy [130]
## $ participant   <chr> "BTCQUV", "BTCQUV", "BTCQUV", "BTCQUV", "BTCQUV", "BTCQU…
## $ session       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ group         <chr> "ws", "ws", "ws", "ws", "ws", "ws", "ws", "ws", "ws", "w…
## $ expName       <chr> "ws-ti-task-prolific", "ws-ti-task-prolific", "ws-ti-tas…
## $ date          <chr> "2022-02-20_16h15.19.643", "2022-02-20_16h15.19.643", "2…
## $ corr          <dbl> 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0,…
## $ rt            <dbl> 1.7423, 1.1811, 1.8598, 1.8475, 2.2589, 1.6005, 3.2540, …
## $ Rank1         <dbl> 4, 2, 5, 1, 3, 5, 4, 2, 6, 3, 3, 6, 4, 5, 3, 2, 5, 2, 4,…
## $ Rank2         <dbl> 5, 3, 4, 2, 4, 6, 3, 1, 5, 2, 2, 5, 3, 6, 4, 1, 4, 3, 5,…
## $ stimCategory  <chr> "Faces", "Faces", "Faces", "Faces", "Faces", "Faces", "F…
## $ part          <fct> learning, learning, learning, learning, learning, learni…
## $ feedback      <dbl> 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0,…
## $ hierarchy     <chr> "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H1", "H…
## $ trial         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
## $ expVersion    <chr> "face-scene", "face-scene", "face-scene", "face-scene", …
## $ participant_u <chr> "BTCQUV_face-scene", "BTCQUV_face-scene", "BTCQUV_face-s…
## $ cRank         <chr> "45", "23", "45", "12", "34", "56", "34", "12", "56", "2…
## $ pairType      <chr> "premise", "premise", "premise", "premise", "premise", "…

# number of participants
n_distinct(exp1.combined$participant_u)

## [1] 70

# experiment structure
DT::datatable(exp1.combined %>% group_by(participant_u, session, part) %>% count())

exp1.combined.imAvg <- exp1.combined %>% filter(part == "immediate testing") %>% group_by(participant_u, hierarchy) %>% summarise(meanPremisePerformance = mean(corr))

## `summarise()` has grouped output by 'participant_u'. You can override using the
## `.groups` argument.

exp1.inf <- exp1.combined %>% filter(part == "delayed testing", pairType=="inference") %>% left_join(exp1.combined.imAvg)

## Joining, by = c("hierarchy", "participant_u")

grouped_ggwithinstats(
  data             = exp1.combined %>% group_by(participant_u, part, pairType, hierarchy) %>% summarise(meanPerf = mean(corr)) %>% filter(pairType != "anchor", part != "learning") %>% mutate(partXpairType = interaction(part, pairType)),
  x                = hierarchy,
  y                = meanPerf,
  grouping.var     = partXpairType,
  type             = "p"
)

## `summarise()` has grouped output by 'participant_u', 'part', 'pairType'. You
## can override using the `.groups` argument.

# Using a mixed logistic model to analyse main effect of Hierarchy (H1: encoded 27h ago, H2: encoded 3h ago)

m1 <- glmer(corr ~ 1 +
    (1 | participant_u), data = exp1.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2 <- glmer(corr ~ hierarchy +
    (1 | participant_u), data = exp1.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3 <- glmer(corr ~ hierarchy + meanPremisePerformance +
    (1 | participant_u), data = exp1.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4 <- glmer(corr ~ hierarchy * meanPremisePerformance +
    (1 | participant_u), data = exp1.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1,m2,m3,m4)

	corr			corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.29	1.00 – 1.66	0.053	1.39	1.06 – 1.82	0.016	0.05	0.03 – 0.11	<0.001	0.02	0.01 – 0.05	<0.001
hierarchy [H2]				0.86	0.74 – 1.00	0.045	0.55	0.46 – 0.65	<0.001	4.65	1.97 – 10.96	<0.001
meanPremisePerformance							101.13	40.23 – 254.26	<0.001	315.68	112.33 – 887.16	<0.001
hierarchy [H2] × meanPremisePerformance										0.06	0.02 – 0.18	<0.001
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.08 _{participant_u}			1.08 _{participant_u}			1.45 _{participant_u}			1.44 _{participant_u}
ICC	0.25			0.25			0.31			0.31
N	70 _{participant_u}			70 _{participant_u}			70 _{participant_u}			70 _{participant_u}
Observations	3360			3360			3360			3360
Marginal R² / Conditional R²	0.000 / 0.247			0.001 / 0.249			0.100 / 0.375			0.105 / 0.378

anova(m1,m2,m3,m4) #m4 wins

## Data: exp1.inf
## Models:
## m1: corr ~ 1 + (1 | participant_u)
## m2: corr ~ hierarchy + (1 | participant_u)
## m3: corr ~ hierarchy + meanPremisePerformance + (1 | participant_u)
## m4: corr ~ hierarchy * meanPremisePerformance + (1 | participant_u)
##    npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
## m1    2 4219.9 4232.1 -2107.9   4215.9                           
## m2    3 4217.9 4236.3 -2105.9   4211.9   4.0046  1    0.04538 *  
## m3    4 4112.5 4137.0 -2052.2   4104.5 107.3867  1  < 2.2e-16 ***
## m4    5 4089.6 4120.2 -2039.8   4079.6  24.9266  1  5.956e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DT::datatable(anova(m1,m2,m3,m4))

plot_model(m4, type = "pred", terms = c("meanPremisePerformance[all]","hierarchy"), axis.title="Inference performance%" )

plot_model(m4, type = "diag")

## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.5.3
## Current Matrix version is 1.5.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package

## $participant_u

## `geom_smooth()` using formula = 'y ~ x'

# Using a mixed logistic model to analyse interaction between hierarchy and distance

exp1.inf_dist <- exp1.combined %>% filter(part == "delayed testing", pairType=="inference") %>% left_join(exp1.combined.imAvg) %>% mutate(distance = abs(Rank1-Rank2), jointrank = Rank1+Rank2)

## Joining, by = c("hierarchy", "participant_u")

#create scatterplot
plot(exp1.inf_dist$distance, exp1.inf_dist$jointrank, main= "Joint rank vs Distance")
#add labels to every point
text(exp1.inf_dist$distance+0.02,exp1.inf_dist$jointrank, labels=exp1.inf_dist$cRank)

m1.d <- glmer(corr ~ hierarchy * meanPremisePerformance +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.d <- glmer(corr ~ hierarchy * meanPremisePerformance + distance +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.d <- glmer(corr ~ hierarchy * meanPremisePerformance + distance * hierarchy +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.d <- glmer(corr ~ hierarchy * meanPremisePerformance * distance +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d,m2.d,m3.d,m4.d)

	corr			corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	0.02	0.01 – 0.05	<0.001	0.01	0.01 – 0.03	<0.001	0.01	0.00 – 0.03	<0.001	0.03	0.00 – 0.28	0.003
hierarchy [H2]	4.65	1.97 – 10.96	<0.001	4.67	1.98 – 11.02	<0.001	6.15	1.97 – 19.24	0.002	16.26	0.40 – 663.23	0.140
meanPremisePerformance	315.68	112.33 – 887.16	<0.001	321.91	114.47 – 905.28	<0.001	322.49	114.64 – 907.21	<0.001	106.33	3.75 – 3012.65	0.006
hierarchy [H2] × meanPremisePerformance	0.06	0.02 – 0.18	<0.001	0.06	0.02 – 0.18	<0.001	0.06	0.02 – 0.18	<0.001	0.02	0.00 – 2.47	0.112
distance				1.29	1.10 – 1.51	0.002	1.37	1.09 – 1.72	0.007	0.98	0.37 – 2.62	0.966
hierarchy [H2] × distance							0.89	0.64 – 1.23	0.472	0.59	0.12 – 2.76	0.500
meanPremisePerformance × distance										1.62	0.41 – 6.37	0.493
(hierarchy [H2] × meanPremisePerformance) × distance										1.59	0.21 – 11.92	0.651
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.44 _{participant_u}			1.45 _{participant_u}			1.45 _{participant_u}			1.46 _{participant_u}
ICC	0.31			0.31			0.31			0.31
N	70 _{participant_u}			70 _{participant_u}			70 _{participant_u}			70 _{participant_u}
Observations	3360			3360			3360			3360
Marginal R² / Conditional R²	0.105 / 0.378			0.108 / 0.381			0.108 / 0.381			0.109 / 0.382

anova(m1.d,m2.d,m3.d,m4.d) # m2 wins

## Data: exp1.inf_dist
## Models:
## m1.d: corr ~ hierarchy * meanPremisePerformance + (1 | participant_u)
## m2.d: corr ~ hierarchy * meanPremisePerformance + distance + (1 | participant_u)
## m3.d: corr ~ hierarchy * meanPremisePerformance + distance * hierarchy + (1 | participant_u)
## m4.d: corr ~ hierarchy * meanPremisePerformance * distance + (1 | participant_u)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1.d    5 4089.6 4120.2 -2039.8   4079.6                        
## m2.d    6 4082.0 4118.7 -2035.0   4070.0 9.5509  1   0.001999 **
## m3.d    7 4083.5 4126.4 -2034.8   4069.5 0.5143  1   0.473288   
## m4.d    9 4085.6 4140.6 -2033.8   4067.6 1.9593  2   0.375445   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m2.d, type = "pred", terms = c("meanPremisePerformance[all]","hierarchy","distance"), axis.title="Inference performance%" )

plot_model(m2.d, type = "diag")

## $participant_u

## `geom_smooth()` using formula = 'y ~ x'

# Using a mixed logistic model to analyse interaction between hierarchy and distance [DISTANT pair only]

m1.d3 <- glmer(corr ~ hierarchy * meanPremisePerformance +
    (1 | participant_u), data = exp1.inf_dist %>% filter(cRank == "25"), family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d3)

	corr
Predictors	Odds Ratios	CI	p
(Intercept)	0.02	0.01 – 0.09	<0.001
hierarchy [H2]	6.28	1.33 – 29.66	0.020
meanPremisePerformance	445.90	72.45 – 2744.29	<0.001
hierarchy [H2] × meanPremisePerformance	0.03	0.00 – 0.26	0.001
Random Effects
σ²	3.29
τ₀₀ _{participant_u}	2.40
ICC	0.42
N _{participant_u}	70
Observations	1120
Marginal R² / Conditional R²	0.097 / 0.478

anova(m1.d3)

## Analysis of Variance Table
##                                  npar Sum Sq Mean Sq F value
## hierarchy                           1  2.770   2.770  2.7703
## meanPremisePerformance              1 36.546  36.546 36.5456
## hierarchy:meanPremisePerformance    1 11.017  11.017 11.0169

plot_model(m2.d, type = "pred", terms = c("meanPremisePerformance[all]","distance"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

m1.j <- glmer(corr ~ hierarchy * meanPremisePerformance +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.j <- glmer(corr ~ hierarchy * meanPremisePerformance + jointrank +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.j <- glmer(corr ~ hierarchy * meanPremisePerformance + jointrank * hierarchy +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.j <- glmer(corr ~ hierarchy * meanPremisePerformance * jointrank +
    (1 | participant_u), data = exp1.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.j,m2.j,m3.j,m4.j)

	corr			corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	0.02	0.01 – 0.05	<0.001	0.03	0.01 – 0.08	<0.001	0.05	0.02 – 0.18	<0.001	0.05	0.00 – 3.07	0.157
hierarchy [H2]	4.65	1.97 – 10.96	<0.001	4.65	1.97 – 10.96	<0.001	1.45	0.31 – 6.86	0.637	5526.47	9.94 – 3073432.49	0.008
meanPremisePerformance	315.68	112.33 – 887.16	<0.001	315.92	112.43 – 887.70	<0.001	318.04	113.09 – 894.40	<0.001	313.78	1.16 – 85123.18	0.044
hierarchy [H2] × meanPremisePerformance	0.06	0.02 – 0.18	<0.001	0.06	0.02 – 0.18	<0.001	0.06	0.02 – 0.18	<0.001	0.00	0.00 – 0.01	0.002
jointrank				0.97	0.89 – 1.07	0.537	0.89	0.78 – 1.02	0.092	0.89	0.50 – 1.57	0.687
hierarchy [H2] × jointrank							1.18	0.98 – 1.42	0.078	0.36	0.15 – 0.89	0.027
meanPremisePerformance × jointrank										1.00	0.46 – 2.21	0.991
(hierarchy [H2] × meanPremisePerformance) × jointrank										4.39	1.38 – 14.02	0.012
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.44 _{participant_u}			1.44 _{participant_u}			1.45 _{participant_u}			1.46 _{participant_u}
ICC	0.31			0.31			0.31			0.31
N	70 _{participant_u}			70 _{participant_u}			70 _{participant_u}			70 _{participant_u}
Observations	3360			3360			3360			3360
Marginal R² / Conditional R²	0.105 / 0.378			0.105 / 0.378			0.106 / 0.379			0.109 / 0.383

anova(m1.j,m2.j,m3.j,m4.j) #m4 wins

## Data: exp1.inf_dist
## Models:
## m1.j: corr ~ hierarchy * meanPremisePerformance + (1 | participant_u)
## m2.j: corr ~ hierarchy * meanPremisePerformance + jointrank + (1 | participant_u)
## m3.j: corr ~ hierarchy * meanPremisePerformance + jointrank * hierarchy + (1 | participant_u)
## m4.j: corr ~ hierarchy * meanPremisePerformance * jointrank + (1 | participant_u)
##      npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)   
## m1.j    5 4089.6 4120.2 -2039.8   4079.6                         
## m2.j    6 4091.2 4127.9 -2039.6   4079.2  0.3786  1   0.538331   
## m3.j    7 4090.1 4132.9 -2038.1   4076.1  3.0885  1   0.078847 . 
## m4.j    9 4083.0 4138.1 -2032.5   4065.0 11.0643  2   0.003958 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compare_performance(m4.d, m4.j)

## # Comparison of Model Performance Indices
## 
## Name |    Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## m4.d | glmerMod | 4085.6 (0.222) | 4085.6 (0.222) | 4140.6 (0.222) |      0.382 |      0.109 | 0.307 | 0.442 | 1.000 |    0.568 |      -Inf |       5.305e-04
## m4.j | glmerMod | 4083.0 (0.778) | 4083.1 (0.778) | 4138.1 (0.778) |      0.383 |      0.109 | 0.307 | 0.442 | 1.000 |    0.568 |      -Inf |       8.418e-04

test_performance(m4.d, m4.j)

## Name |    Model |   BF
## ----------------------
## m4.d | glmerMod |     
## m4.j | glmerMod | 3.50
## Each model is compared to m4.d.

plot(compare_performance(m4.d, m4.j))

plot_model(m4.j, type = "pred", terms = c("jointrank","hierarchy"), axis.title="Inference performance%" )

plot_model(m4.j, type = "pred", terms = c("meanPremisePerformance[all]","jointrank", "hierarchy"), axis.title="Inference performance%" )

plot_model(m4.j, type = "diag")

## $participant_u

## `geom_smooth()` using formula = 'y ~ x'

grouped_ggwithinstats(
  data             = exp1.inf_dist %>% group_by(participant_u, cRank, hierarchy) %>% summarise(meanPerf = mean(corr)),
  x                = hierarchy,
  y                = meanPerf,
  grouping.var     = cRank,
  type             = "p"
)

## `summarise()` has grouped output by 'participant_u', 'cRank'. You can override
## using the `.groups` argument.

Dataset 2 (Between-subject Transitive Inference with 12h delayed testing) by Lorena n=24 (only condition PRE)

best fitting model is the one that includes interaction between Group ( time of encoding ) and Immediate premise pair performance ( strength of memory encoding ). Interaction is significant.
main effect of distance, and significant distance and time of encoding interaction!
- but a joint rank of items model fits the data sig better.
there is a 3way interaction between strength of memory encoding, time of encoding and joint rank of items hinting at a certain kind of “spatial” time-dependent effect on inference performance
- the time of encoding and joint rank of items two-way interaction is also sig.

ellrep.lab.learning <- read_csv("lorena-ti-pilot-ellrep/Learning_Beh2.csv") %>% filter(Condition=="pre") %>% mutate(part = "learning", RT = as.character(RT))

## Rows: 8730 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (7): ID, Group, Condition, Pair, Stim, session, PairType
## dbl (4): acc, RT, Distance, Block
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ellrep.lab.im <- read_delim("lorena-ti-pilot-ellrep/T_premises_pre - copia.csv", 
    delim = ";", escape_double = FALSE, trim_ws = TRUE) %>% group_by(ID_S1,Group_S1,Stim_S1)

## Rows: 2880 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (16): ID_S1, Group_S1, Condition_S1, RT_S1, Pair_S1, Stim_S1, session_S1...
## dbl  (6): acc_S1, Distance_S1, Block_S1, acc_S2, Distance_S2, Block_S2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ellrep.lab.im.aggr <- read_delim("lorena-ti-pilot-ellrep/T_premises_pre - copia.csv", 
    delim = ";", escape_double = FALSE, trim_ws = TRUE) %>% group_by(ID_S1,Group_S1,Stim_S1) %>% summarise(meanIm = mean(acc_S1))

## Rows: 2880 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (16): ID_S1, Group_S1, Condition_S1, RT_S1, Pair_S1, Stim_S1, session_S1...
## dbl  (6): acc_S1, Distance_S1, Block_S1, acc_S2, Distance_S2, Block_S2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## `summarise()` has grouped output by 'ID_S1', 'Group_S1'. You can override using the `.groups` argument.

read_csv("lorena-ti-pilot-ellrep/Learning_Beh2.csv") %>% group_by(Group,Condition) %>% summarise(n_distinct(ID))

## Rows: 8730 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (7): ID, Group, Condition, Pair, Stim, session, PairType
## dbl (4): acc, RT, Distance, Block
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## `summarise()` has grouped output by 'Group'. You can override using the `.groups` argument.

## # A tibble: 4 × 3
## # Groups:   Group [2]
##   Group Condition `n_distinct(ID)`
##   <chr> <chr>                <int>
## 1 sleep all                     11
## 2 sleep pre                     12
## 3 Wake  all                     11
## 4 Wake  pre                     12

ellrep.lab <- read_delim("lorena-ti-pilot-ellrep/T_premises_pre - copia.csv",delim =';') %>% pivot_longer(
    cols = contains("_"), 
    names_to = c('.value', 'session'),
    values_drop_na = TRUE,
    names_pattern = '(.*)\\_(S\\d+)',
    names_repair = "unique"
  ) %>% rename(session = "session...1") %>% select(-session...9) -> ellrep.lab.tmp

## Rows: 2880 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (16): ID_S1, Group_S1, Condition_S1, RT_S1, Pair_S1, Stim_S1, session_S1...
## dbl  (6): acc_S1, Distance_S1, Block_S1, acc_S2, Distance_S2, Block_S2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## New names:

ellrep.lab <- rbind(ellrep.lab ,read_delim("lorena-ti-pilot-ellrep/T_inferences_pre - copia.csv",delim =';'),read_delim("lorena-ti-pilot-ellrep/T_anchor_pre - copia.csv",delim =';'))

## Rows: 1728 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (8): ID, Group, Condition, RT, Pair, Stim, session, PairType
## dbl (3): acc, Distance, Block
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 576 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (8): ID, Group, Condition, RT, Pair, Stim, session, PairType
## dbl (3): acc, Distance, Block
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ellrep.lab <- bind_rows(ellrep.lab,ellrep.lab.learning) %>% rowwise() %>% mutate(part = ifelse(session == "S1" && is.na(part), "immediate testing",part)) %>% mutate(part = ifelse(session == "S2" && is.na(part), "delayed testing",part))

# Split the "pair" column into two separate columns
split_df <- lapply(strsplit(ellrep.lab$Pair, ""), function(x) x[1:2])
ellrep.lab[c("Rank1", "Rank2")] <- data.frame(do.call(rbind, split_df))
# Convert each letter to its corresponding position in the alphabet
ellrep.lab$Rank1 <- as.numeric(chartr("ABCDEFGHIJKLMNOPQRSTUVWXYZ", "1234567891011121314151617181920212223242526", ellrep.lab$Rank1))
ellrep.lab$Rank2 <- as.numeric(chartr("ABCDEFGHIJKLMNOPQRSTUVWXYZ", "1234567891011121314151617181920212223242526", ellrep.lab$Rank2))

ellrep.lab %>% rename(participant=ID, group=Group, pair=Pair, pairType=PairType,stimCategory=Stim) %>% mutate(session = factor( readr::parse_number(as.character(session))),stimCategory = ifelse(stimCategory=="F","faces",ifelse(stimCategory=="O","objects","scenes")), pairType=tolower(pairType), group=tolower(group), participant=tolower(participant), distance = abs(Rank1-Rank2), jointrank=Rank1+Rank2) -> exp2.combined

glimpse(exp2.combined)

## Rows: 12,384
## Columns: 16
## Rowwise: 
## $ session      <fct> 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, …
## $ participant  <chr> "participant_1_sleep", "participant_1_sleep", "participan…
## $ group        <chr> "sleep", "sleep", "sleep", "sleep", "sleep", "sleep", "sl…
## $ Condition    <chr> "pre", "pre", "pre", "pre", "pre", "pre", "pre", "pre", "…
## $ acc          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
## $ RT           <chr> "139.316.063.878.141", "171.018.139.833.234", "170.355.61…
## $ pair         <chr> "AB", "AB", "AB", "AB", "BC", "BC", "BC", "BC", "CD", "CD…
## $ stimCategory <chr> "faces", "faces", "faces", "faces", "faces", "faces", "fa…
## $ pairType     <chr> "premise", "premise", "premise", "premise", "premise", "p…
## $ Distance     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Block        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ part         <chr> "immediate testing", "delayed testing", "immediate testin…
## $ Rank1        <dbl> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, …
## $ Rank2        <dbl> 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, …
## $ distance     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ jointrank    <dbl> 3, 3, 3, 3, 5, 5, 5, 5, 7, 7, 7, 7, 9, 9, 9, 9, 11, 11, 1…

exp2.combined %>% group_by(participant) %>% count()

## # A tibble: 27 × 2
## # Groups:   participant [27]
##    participant              n
##    <chr>                <int>
##  1 participant_1_sleep    516
##  2 participant_1_wake     606
##  3 participant_11_wake    456
##  4 participant_12_wake    636
##  5 participant_13_sleep   426
##  6 participant_13_wake    486
##  7 participant_14_sleep   516
##  8 participant_14_wake    606
##  9 participant_15_sleep   456
## 10 participant_15_wake    576
## # … with 17 more rows

#write.csv(exp2.combined,"ellrep-lab_bylorena.csv")

# number of participants
n_distinct(exp2.combined$participant)

## [1] 27

# experiment structure
DT::datatable(exp2.combined %>% group_by(group, participant, session, part, stimCategory) %>% count())

grouped_ggbetweenstats(
  data             = exp2.combined %>% group_by(participant, group, pairType, part) %>% summarise(meanPerf = mean(acc)) %>% filter(pairType != "anchor", part != "learning") %>% mutate(partXpairType = interaction(part, pairType))  %>% mutate(partXpairType = factor(interaction(part, pairType), levels=c("immediate testing.premise", "delayed testing.premise", "delayed testing.inference"))),
  x                = group,
  y                = meanPerf,
  grouping.var     = partXpairType,
  type             = "p",
    bf.message = FALSE
)

## `summarise()` has grouped output by 'participant', 'group', 'pairType'. You can
## override using the `.groups` argument.

grouped_ggbetweenstats(
  data             = exp2.combined %>% group_by(participant, group, pairType, part, stimCategory) %>% summarise(meanPerf = mean(acc)) %>% filter(pairType != "anchor", part != "learning") %>% mutate(partXpairType = factor(interaction(part, pairType), levels=c("immediate testing.premise", "delayed testing.premise", "delayed testing.inference"))),
  x                = group,
  y                = meanPerf,
  grouping.var     = partXpairType,
  type             = "p",
  bf.message = FALSE
)

## `summarise()` has grouped output by 'participant', 'group', 'pairType', 'part'.
## You can override using the `.groups` argument.

exp2.combined.imAvg <- exp2.combined %>% filter(part == "immediate testing") %>% group_by(participant, group, stimCategory) %>% summarise(meanPremisePerformance = mean(acc))

## `summarise()` has grouped output by 'participant', 'group'. You can override
## using the `.groups` argument.

exp2.inf <- exp2.combined %>% filter(part == "delayed testing", pairType=="inference") %>% left_join(exp2.combined.imAvg)

## Joining, by = c("participant", "group", "stimCategory")

No random slope

# Using a mixed logistic model to analyse main effect of ...

m1 <- glmer(acc ~ 1 +
    (1 | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2 <- glmer(acc ~ group +
    (1 | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3 <- glmer(acc ~ group + meanPremisePerformance +
    (1 | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4 <- glmer(acc ~ group * meanPremisePerformance +
    (1 | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1,m2,m3,m4)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.53	1.00 – 2.34	0.049	2.04	1.14 – 3.67	0.017	3.99	1.75 – 9.07	0.001	1.05	0.38 – 2.91	0.932
group [wake]				0.56	0.25 – 1.28	0.170	0.57	0.24 – 1.35	0.199	4.98	1.28 – 19.30	0.020
meanPremisePerformance							0.38	0.18 – 0.83	0.015	2.64	0.77 – 8.99	0.121
group [wake] × meanPremisePerformance										0.04	0.01 – 0.21	<0.001
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.05 _participant			0.97 _participant			1.10 _participant			0.93 _participant
ICC	0.24			0.23			0.25			0.22
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.000 / 0.241			0.019 / 0.242			0.028 / 0.271			0.050 / 0.258

plot_model(m4, type = "pred", terms = c("meanPremisePerformance[all]","group"), axis.title="Inference performance%" )

# Using a mixed logistic model to analyse interaction between hierarchy and distance

exp2.inf_dist <- exp2.combined %>% filter(part == "delayed testing", pairType=="inference") %>% left_join(exp2.combined.imAvg) %>% mutate(distance = abs(Rank1-Rank2), jointrank = Rank1+Rank2, cRank = paste0(Rank1,Rank2))

## Joining, by = c("participant", "group", "stimCategory")

m1.d <- glmer(acc ~ group * meanPremisePerformance +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.d <- glmer(acc ~ group * meanPremisePerformance + distance +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.d <- glmer(acc ~ group * meanPremisePerformance + distance * group +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.d <- glmer(acc ~ group * meanPremisePerformance * distance +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d,m2.d,m3.d,m4.d)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.05	0.38 – 2.91	0.932	0.96	0.30 – 3.01	0.941	0.47	0.13 – 1.69	0.249	0.12	0.01 – 1.95	0.137
group [wake]	4.98	1.28 – 19.30	0.020	4.98	1.28 – 19.29	0.020	18.94	3.42 – 104.79	0.001	11.72	0.30 – 462.62	0.189
meanPremisePerformance	2.64	0.77 – 8.99	0.121	2.64	0.77 – 8.99	0.121	2.65	0.77 – 9.06	0.120	20.57	0.41 – 1037.19	0.131
group [wake] × meanPremisePerformance	0.04	0.01 – 0.21	<0.001	0.04	0.01 – 0.21	<0.001	0.04	0.01 – 0.21	<0.001	0.08	0.00 – 12.44	0.324
distance				1.04	0.83 – 1.29	0.736	1.41	1.02 – 1.94	0.039	2.51	0.83 – 7.54	0.102
group [wake] × distance							0.57	0.36 – 0.88	0.011	0.69	0.16 – 3.02	0.627
meanPremisePerformance × distance										0.41	0.08 – 2.05	0.280
(group [wake] × meanPremisePerformance) × distance										0.78	0.10 – 6.23	0.816
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	0.93 _participant			0.93 _participant			0.93 _participant			0.93 _participant
ICC	0.22			0.22			0.22			0.22
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.050 / 0.258			0.050 / 0.258			0.054 / 0.263			0.056 / 0.265

anova(m1.d,m2.d,m3.d,m4.d) #m3 wins

## Data: exp2.inf_dist
## Models:
## m1.d: acc ~ group * meanPremisePerformance + (1 | participant)
## m2.d: acc ~ group * meanPremisePerformance + distance + (1 | participant)
## m3.d: acc ~ group * meanPremisePerformance + distance * group + (1 | participant)
## m4.d: acc ~ group * meanPremisePerformance * distance + (1 | participant)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1.d    5 2136.3 2163.5 -1063.1   2126.3                       
## m2.d    6 2138.1 2170.9 -1063.1   2126.1 0.1132  1    0.73657  
## m3.d    7 2133.7 2171.9 -1059.9   2119.7 6.4085  1    0.01136 *
## m4.d    9 2133.8 2182.9 -1057.9   2115.8 3.9345  2    0.13984  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m3.d, type = "pred", terms = c("group","distance"), axis.title="Inference performance%" )

plot_model(m3.d, type = "pred", terms = c("meanPremisePerformance[all]","group","distance"), axis.title="Inference performance%" )

plot_model(m3.d, type = "diag")

## $participant

## `geom_smooth()` using formula = 'y ~ x'

# Using a mixed logistic model to analyse interaction between hierarchy and distance [DISTANT pair only]

m1.d3 <- glmer(acc ~ group * meanPremisePerformance +
    (1 | participant), data = exp2.inf_dist %>% filter(distance == "3"), family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d3)

	acc
Predictors	Odds Ratios	CI	p
(Intercept)	2.37	0.43 – 13.11	0.321
group [wake]	3.02	0.33 – 27.56	0.327
meanPremisePerformance	1.47	0.17 – 12.64	0.727
group [wake] × meanPremisePerformance	0.04	0.00 – 0.61	0.021
Random Effects
σ²	3.29
τ₀₀ _participant	1.61
ICC	0.33
N _participant	24
Observations	576
Marginal R² / Conditional R²	0.102 / 0.397

m1.j <- glmer(acc ~ group * meanPremisePerformance +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.j <- glmer(acc ~ group * meanPremisePerformance + jointrank +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.j <- glmer(acc ~ group * meanPremisePerformance + jointrank * group +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.j <- glmer(acc ~ group * meanPremisePerformance * jointrank +
    (1 | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.j,m2.j,m3.j,m4.j)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.05	0.38 – 2.91	0.932	1.36	0.35 – 5.28	0.655	2.50	0.48 – 13.00	0.277	0.00	0.00 – 0.27	0.012
group [wake]	4.98	1.28 – 19.30	0.020	4.98	1.29 – 19.30	0.020	1.58	0.17 – 14.75	0.689	415.68	0.92 – 188676.08	0.053
meanPremisePerformance	2.64	0.77 – 8.99	0.121	2.64	0.77 – 8.99	0.121	2.64	0.77 – 9.02	0.121	93674.98	108.23 – 81079248.72	0.001
group [wake] × meanPremisePerformance	0.04	0.01 – 0.21	<0.001	0.04	0.01 – 0.21	<0.001	0.04	0.01 – 0.21	<0.001	0.00	0.00 – 0.04	0.008
jointrank				0.96	0.85 – 1.09	0.560	0.88	0.73 – 1.06	0.187	2.35	1.23 – 4.47	0.009
group [wake] × jointrank							1.18	0.91 – 1.52	0.205	0.53	0.23 – 1.25	0.147
meanPremisePerformance × jointrank										0.22	0.09 – 0.58	0.002
(group [wake] × meanPremisePerformance) × jointrank										3.42	1.02 – 11.48	0.046
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	0.93 _participant			0.93 _participant			0.93 _participant			0.96 _participant
ICC	0.22			0.22			0.22			0.23
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.050 / 0.258			0.050 / 0.259			0.051 / 0.260			0.060 / 0.272

anova(m1.j,m2.j,m3.j,m4.j) #m4 wins

## Data: exp2.inf_dist
## Models:
## m1.j: acc ~ group * meanPremisePerformance + (1 | participant)
## m2.j: acc ~ group * meanPremisePerformance + jointrank + (1 | participant)
## m3.j: acc ~ group * meanPremisePerformance + jointrank * group + (1 | participant)
## m4.j: acc ~ group * meanPremisePerformance * jointrank + (1 | participant)
##      npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)   
## m1.j    5 2136.3 2163.5 -1063.1   2126.3                         
## m2.j    6 2137.9 2170.6 -1063.0   2125.9  0.3394  1   0.560174   
## m3.j    7 2138.3 2176.5 -1062.2   2124.3  1.5995  1   0.205971   
## m4.j    9 2132.1 2181.2 -1057.1   2114.1 10.2001  2   0.006096 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m4.j, type = "pred", terms = c("jointrank","group"), axis.title="Inference performance%" )

plot_model(m4.j, type = "pred", terms = c("meanPremisePerformance[all]","jointrank", "group"), axis.title="Inference performance%" )

plot_model(m4.j, type = "diag")

## $participant

## `geom_smooth()` using formula = 'y ~ x'

# distance model vs jointrank model? jointrank model outperforms...
anova(m3.d,m4.j)

## Data: exp2.inf_dist
## Models:
## m3.d: acc ~ group * meanPremisePerformance + distance * group + (1 | participant)
## m4.j: acc ~ group * meanPremisePerformance * jointrank + (1 | participant)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m3.d    7 2133.7 2171.9 -1059.9   2119.7                       
## m4.j    9 2132.1 2181.2 -1057.1   2114.1 5.6173  2    0.06029 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compare_performance(m3.d,m4.d, m4.j)

## # Comparison of Model Performance Indices
## 
## Name |    Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## m3.d | glmerMod | 2133.7 (0.237) | 2133.8 (0.241) | 2171.9 (0.986) |      0.263 |      0.054 | 0.221 | 0.452 | 1.000 |    0.589 |      -Inf |           0.002
## m4.d | glmerMod | 2133.8 (0.230) | 2133.9 (0.229) | 2182.9 (0.004) |      0.265 |      0.056 | 0.221 | 0.452 | 1.000 |    0.588 |      -Inf |           0.002
## m4.j | glmerMod | 2132.1 (0.533) | 2132.2 (0.530) | 2181.2 (0.009) |      0.272 |      0.060 | 0.225 | 0.451 | 1.000 |    0.587 |      -Inf |           0.002

plot(compare_performance(m3.d,m4.d, m4.j))

test_performance(m3.d,m4.d, m4.j)

## Name |    Model |    BF
## -----------------------
## m3.d | glmerMod |      
## m4.d | glmerMod | 0.004
## m4.j | glmerMod | 0.010
## Each model is compared to m3.d.

With random slope

m1 <- glmer(acc ~ 1 +
    (stimCategory | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2 <- glmer(acc ~ group +
    (stimCategory| participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3 <- glmer(acc ~ group + meanPremisePerformance +
    (stimCategory | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4 <- glmer(acc ~ group * meanPremisePerformance +
    (stimCategory | participant), data = exp2.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1,m2,m3,m4)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.70	0.99 – 2.94	0.056	2.56	1.25 – 5.27	0.011	2.78	0.41 – 18.70	0.293	0.25	0.02 – 2.86	0.262
group [wake]				0.43	0.15 – 1.18	0.101	0.43	0.15 – 1.18	0.101	14.42	0.82 – 253.54	0.068
meanPremisePerformance							0.89	0.07 – 11.28	0.928	29.07	0.97 – 874.22	0.052
group [wake] × meanPremisePerformance										0.01	0.00 – 0.31	0.011
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	2.96 _participant			2.76 _participant			2.80 _participant			1.87 _participant
τ₁₁	4.39 _{participant.stimCategoryobjects}			4.46 _{participant.stimCategoryobjects}			4.42 _{participant.stimCategoryobjects}			4.97 _{participant.stimCategoryobjects}
	4.83 _{participant.stimCategoryscenes}			4.82 _{participant.stimCategoryscenes}			4.81 _{participant.stimCategoryscenes}			4.29 _{participant.stimCategoryscenes}
ρ₀₁	-0.51			-0.47			-0.47			-0.47
	-0.62			-0.66			-0.66			-0.59
ICC	0.47			0.46			0.46			0.36
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.000 / 0.474			0.029 / 0.472			0.029 / 0.475			0.089 / 0.420

anova(m1,m2,m3,m4) #m4 wins

## Data: exp2.inf
## Models:
## m1: acc ~ 1 + (stimCategory | participant)
## m2: acc ~ group + (stimCategory | participant)
## m3: acc ~ group + meanPremisePerformance + (stimCategory | participant)
## m4: acc ~ group * meanPremisePerformance + (stimCategory | participant)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1    7 1898.5 1936.6 -942.23   1884.5                       
## m2    8 1897.9 1941.5 -940.95   1881.9 2.5680  1     0.1090  
## m3    9 1899.9 1949.0 -940.94   1881.9 0.0079  1     0.9291  
## m4   10 1896.2 1950.8 -938.10   1876.2 5.6761  1     0.0172 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m4, type = "pred", terms = c("meanPremisePerformance[all]","group"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m4, type = "diag")

## $participant

## `geom_smooth()` using formula = 'y ~ x'

# Using a mixed logistic model to analyse interaction between hierarchy and distance

exp2.inf_dist <- exp2.combined %>% filter(part == "delayed testing", pairType=="inference") %>% left_join(exp2.combined.imAvg) %>% mutate(distance = abs(Rank1-Rank2), jointrank = Rank1+Rank2, cRank = paste0(Rank1,Rank2))

## Joining, by = c("participant", "group", "stimCategory")

m1.d <- glmer(acc ~ group * meanPremisePerformance +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.d <- glmer(acc ~ group * meanPremisePerformance + distance +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.d <- glmer(acc ~ group * meanPremisePerformance + distance * group +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.d <- glmer(acc ~ group * meanPremisePerformance * distance +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d,m2.d,m3.d,m4.d)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	0.25	0.02 – 2.86	0.262	0.22	0.02 – 2.78	0.242	0.09	0.01 – 1.26	0.074	0.02	0.00 – 0.83	0.040
group [wake]	14.42	0.82 – 253.54	0.068	14.42	0.80 – 259.60	0.070	77.57	3.38 – 1779.89	0.006	33.94	0.31 – 3710.46	0.141
meanPremisePerformance	29.07	0.97 – 874.22	0.052	29.08	0.94 – 899.38	0.054	29.77	0.94 – 940.77	0.054	307.75	1.50 – 62950.02	0.035
group [wake] × meanPremisePerformance	0.01	0.00 – 0.31	0.011	0.01	0.00 – 0.32	0.012	0.01	0.00 – 0.32	0.012	0.02	0.00 – 13.49	0.236
distance				1.05	0.82 – 1.34	0.703	1.53	1.07 – 2.18	0.020	2.94	0.89 – 9.71	0.078
group [wake] × distance							0.49	0.30 – 0.80	0.004	0.70	0.14 – 3.41	0.654
meanPremisePerformance × distance										0.37	0.06 – 2.09	0.258
(group [wake] × meanPremisePerformance) × distance										0.62	0.06 – 6.01	0.680
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.87 _participant			1.87 _participant			1.89 _participant			1.92 _participant
τ₁₁	4.97 _{participant.stimCategoryobjects}			4.97 _{participant.stimCategoryobjects}			5.03 _{participant.stimCategoryobjects}			5.06 _{participant.stimCategoryobjects}
	4.29 _{participant.stimCategoryscenes}			4.29 _{participant.stimCategoryscenes}			4.34 _{participant.stimCategoryscenes}			4.39 _{participant.stimCategoryscenes}
ρ₀₁	-0.47			-0.47			-0.47			-0.47
	-0.59			-0.59			-0.59			-0.59
ICC	0.36			0.36			0.37			0.37
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.089 / 0.420			0.089 / 0.420			0.095 / 0.425			0.097 / 0.429

anova(m1.d,m2.d,m3.d,m4.d) #m3 wins

## Data: exp2.inf_dist
## Models:
## m1.d: acc ~ group * meanPremisePerformance + (stimCategory | participant)
## m2.d: acc ~ group * meanPremisePerformance + distance + (stimCategory | participant)
## m3.d: acc ~ group * meanPremisePerformance + distance * group + (stimCategory | participant)
## m4.d: acc ~ group * meanPremisePerformance * distance + (stimCategory | participant)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1.d   10 1896.2 1950.8 -938.10   1876.2                        
## m2.d   11 1898.1 1958.1 -938.03   1876.1 0.1421  1   0.706193   
## m3.d   12 1892.0 1957.5 -934.02   1868.0 8.0149  1   0.004639 **
## m4.d   14 1890.9 1967.3 -931.45   1862.9 5.1525  2   0.076058 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m3.d, type = "pred", terms = c("group","distance"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m3.d, type = "pred", terms = c("meanPremisePerformance[all]","group","distance"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m3.d, type = "diag")

## $participant

## `geom_smooth()` using formula = 'y ~ x'

# Using a mixed logistic model to analyse interaction between hierarchy and distance [DISTANT pair only]

m1.d3 <- glmer(acc ~ group * meanPremisePerformance +
    (stimCategory | participant), data = exp2.inf_dist %>% filter(distance == "3"), family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1.d3)

	acc
Predictors	Odds Ratios	CI	p
(Intercept)	0.42	0.01 – 17.25	0.647
group [wake]	29.34	0.25 – 3413.09	0.164
meanPremisePerformance	32.97	0.21 – 5166.14	0.175
group [wake] × meanPremisePerformance	0.00	0.00 – 0.55	0.032
Random Effects
σ²	3.29
τ₀₀ _participant	3.89
τ₁₁ _{participant.stimCategoryobjects}	9.59
τ₁₁ _{participant.stimCategoryscenes}	9.11
ρ₀₁	-0.51
	-0.59
ICC	0.54
N _participant	24
Observations	576
Marginal R² / Conditional R²	0.143 / 0.607

m1.j <- glmer(acc ~ group * meanPremisePerformance +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.j <- glmer(acc ~ group * meanPremisePerformance + jointrank +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3.j <- glmer(acc ~ group * meanPremisePerformance + jointrank * group +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4.j <- glmer(acc ~ group * meanPremisePerformance * jointrank +
    (stimCategory | participant), data = exp2.inf_dist, family = binomial, control = glmerControl(optimizer = "bobyqa"))

## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations
## exceeded

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00515807 (tol = 0.002, component 1)

tab_model(m1.j,m2.j,m3.j,m4.j)

	acc			acc			acc			acc
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	0.25	0.02 – 2.86	0.262	0.34	0.02 – 4.87	0.429	0.72	0.04 – 12.39	0.822	0.00	0.00 – 0.04	0.002
group [wake]	14.42	0.82 – 253.54	0.068	14.43	0.81 – 257.45	0.069	3.46	0.11 – 112.76	0.485	4010.18	3.01 – 5337667.37	0.024
meanPremisePerformance	29.07	0.97 – 874.22	0.052	29.10	0.95 – 892.63	0.054	29.34	0.95 – 904.44	0.053	16244100.05	3964.85 – 66552456734.09	<0.001
group [wake] × meanPremisePerformance	0.01	0.00 – 0.31	0.011	0.01	0.00 – 0.32	0.011	0.01	0.00 – 0.31	0.011	0.00	0.00 – 0.00	0.002
jointrank				0.95	0.83 – 1.10	0.509	0.86	0.70 – 1.05	0.137	2.87	1.41 – 5.85	0.004
group [wake] × jointrank							1.23	0.93 – 1.63	0.154	0.45	0.18 – 1.15	0.095
meanPremisePerformance × jointrank										0.15	0.05 – 0.44	0.001
(group [wake] × meanPremisePerformance) × jointrank										4.77	1.24 – 18.34	0.023
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	1.87 _participant			1.87 _participant			1.88 _participant			1.89 _participant
τ₁₁	4.97 _{participant.stimCategoryobjects}			4.97 _{participant.stimCategoryobjects}			4.99 _{participant.stimCategoryobjects}			5.12 _{participant.stimCategoryobjects}
	4.29 _{participant.stimCategoryscenes}			4.29 _{participant.stimCategoryscenes}			4.30 _{participant.stimCategoryscenes}			4.40 _{participant.stimCategoryscenes}
ρ₀₁	-0.47			-0.47			-0.47			-0.46
	-0.59			-0.59			-0.59			-0.58
ICC	0.36			0.36			0.36			0.36
N	24 _participant			24 _participant			24 _participant			24 _participant
Observations	1728			1728			1728			1728
Marginal R² / Conditional R²	0.089 / 0.420			0.090 / 0.420			0.091 / 0.421			0.106 / 0.432

anova(m1.j,m2.j,m3.j,m4.j) #m4 wins

## Data: exp2.inf_dist
## Models:
## m1.j: acc ~ group * meanPremisePerformance + (stimCategory | participant)
## m2.j: acc ~ group * meanPremisePerformance + jointrank + (stimCategory | participant)
## m3.j: acc ~ group * meanPremisePerformance + jointrank * group + (stimCategory | participant)
## m4.j: acc ~ group * meanPremisePerformance * jointrank + (stimCategory | participant)
##      npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)   
## m1.j   10 1896.2 1950.8 -938.10   1876.2                         
## m2.j   11 1897.8 1957.8 -937.89   1875.8  0.4263  1    0.51383   
## m3.j   12 1897.8 1963.2 -936.90   1873.8  1.9852  1    0.15885   
## m4.j   14 1888.9 1965.2 -930.44   1860.9 12.9139  2    0.00157 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m4.j, type = "pred", terms = c("jointrank","group"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m4.j, type = "pred", terms = c("meanPremisePerformance[all]","jointrank", "group"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m4.j, type = "diag")

## $participant

## `geom_smooth()` using formula = 'y ~ x'

# distance model vs jointrank model? jointrank model outperforms...
anova(m3.d,m4.j)

## Data: exp2.inf_dist
## Models:
## m3.d: acc ~ group * meanPremisePerformance + distance * group + (stimCategory | participant)
## m4.j: acc ~ group * meanPremisePerformance * jointrank + (stimCategory | participant)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m3.d   12 1892.0 1957.5 -934.02   1868.0                       
## m4.j   14 1888.9 1965.2 -930.44   1860.9 7.1683  2    0.02776 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compare_performance(m3.d,m4.d, m4.j)

## # Comparison of Model Performance Indices
## 
## Name |    Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## m3.d | glmerMod | 1892.0 (0.131) | 1892.2 (0.134) | 1957.5 (0.972) |      0.425 |      0.095 | 0.365 | 0.399 | 1.000 |    0.475 |      -Inf |           0.002
## m4.d | glmerMod | 1890.9 (0.232) | 1891.1 (0.231) | 1967.3 (0.007) |      0.429 |      0.097 | 0.368 | 0.398 | 1.000 |    0.474 |      -Inf |           0.002
## m4.j | glmerMod | 1888.9 (0.637) | 1889.1 (0.634) | 1965.2 (0.020) |      0.432 |      0.106 | 0.365 | 0.397 | 1.000 |    0.473 |      -Inf |           0.003

plot(compare_performance(m3.d,m4.d, m4.j))

test_performance(m3.d,m4.d, m4.j)

## Name |    Model |    BF
## -----------------------
## m3.d | glmerMod |      
## m4.d | glmerMod | 0.008
## m4.j | glmerMod | 0.021
## Each model is compared to m3.d.

grouped_ggbetweenstats(
  data             = exp2.inf_dist %>% group_by(participant, pair, group) %>% summarise(meanPerf = mean(acc)),
  x                = group,
  y                = meanPerf,
  grouping.var     = pair,
  type             = "p"
)

## `summarise()` has grouped output by 'participant', 'pair'. You can override
## using the `.groups` argument.

grouped_ggbetweenstats(
  data             = exp2.inf_dist %>% group_by(participant, pair,stimCategory, group) %>% summarise(meanPerf = mean(acc)),
  x                = group,
  y                = meanPerf,
  grouping.var     = pair,
  type             = "p"
)

## `summarise()` has grouped output by 'participant', 'pair', 'stimCategory'. You
## can override using the `.groups` argument.

Dataset 3: (Within-subject Transitive Inference with 24h and 0h delayed testing) - Reanalysis of Berens & Bird (2022) n=17

7 item hierarchy!
Participants trained to ceiling! >> No learning and immediate testing data, trained 4 times longer
main effect of Hierarchy ( time of encoding ) but in the opposite direction: participants were better on inferences on the recently learned hierarchy.
sig distance main effect but in the opposite direction: participants showed less distance effect for the remote hierarchy.
sig joint rank of items and time of encoding interaction but in the opposite direction: recently learned hierarchy show the jointrank effect more and here participants are better for high rank vs low rank items.
Overall distance and joint rank of items similarly good predictors in this dataset, hard to disentangle; as far as time-dependent changes are considered joint rank still seems to show.

exp3 <- read_csv("behrensnbird/DataTable01.csv") %>% mutate(cPairId= paste0(substr(PairId,1,1),parse_number(PairId))) %>% filter(Progressive==0)

## Rows: 816 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): PairId, SubjectId
## dbl (5): k, n, Inferred, Remote, Progressive
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Create a lookup table
lookup_table <- data.frame(
  cPairId = c("D1", "D2", "D3", "D4", "D5", "D6", "I1", "I2", "I3", "I4", "I5", "I6"),
  cPairLetters = c("AB", "BC", "CD", "DE", "EF", "FG", "BD", "BE", "BF", "CE", "CF", "DF")
)
# Use match function to create an index vector
index <- match(exp3$cPairId, lookup_table$cPairId)
# Use index vector to replace values in the column
exp3$cPairLetters <- lookup_table$cPairLetters[index]

# Use str_split_fixed to split the letters column into two separate columns
df_split <- str_split_fixed(exp3$cPairLetters, "", 2)
# Rename the columns in the new dataframe
colnames(df_split) <- c("Rank1", "Rank2")
# Combine the original and new dataframes into a single dataframe
exp3 <- cbind(exp3, df_split)
letter_pairs_to_numbers <- function(x) {
  sapply(strsplit(x, ""), function(y) sum(as.numeric(charToRaw(y)) - 64))
}
exp3$Rank1 <- letter_pairs_to_numbers(exp3$Rank1)
exp3$Rank2 <- letter_pairs_to_numbers(exp3$Rank2)

exp3 <- exp3 %>% mutate(distance = abs(Rank1-Rank2), jointrank = Rank1+Rank2)

# Define a function to convert the k column to binary rows
k_to_binary_rows <- function(x) {
  binary_rows <- matrix(0, nrow = 8, ncol = 1)
  if (x!=0) {
    binary_rows[1:x, 1] <- 1
  }
  return(t(binary_rows))
}

# Apply the function to the k column
df_binary <- as.data.frame(do.call(rbind, lapply(exp3$k, k_to_binary_rows)))

# Rename the columns in the new dataframe
colnames(df_binary) <- paste0("k_", 1:8)

# Combine the original and new dataframes into a single dataframe
exp3_ext <- cbind(exp3, df_binary)

# Use pivot_longer to convert the k columns into a long format
exp3_ext <- pivot_longer(exp3_ext, starts_with("k_"), names_to = "k_ind", values_to = "corr") %>% mutate(pairType = ifelse(Inferred==0,"premise","inference"), hierarchy = ifelse(Remote==1, "H1","H2"))

n_distinct(exp3_ext$SubjectId)

## [1] 17

grouped_ggwithinstats(
  data             = exp3_ext %>% group_by(SubjectId, pairType, hierarchy) %>% summarise(meanPerf = mean(corr)),
  x                = hierarchy,
  y                = meanPerf,
  grouping.var = pairType,
  type             = "p",
  bf.message = FALSE
)

## `summarise()` has grouped output by 'SubjectId', 'pairType'. You can override
## using the `.groups` argument.

exp3.inf <- exp3_ext %>% filter(Inferred==1)

grouped_ggwithinstats(
  data             = exp3.inf %>% group_by(SubjectId, cPairLetters, hierarchy) %>% summarise(meanPerf = mean(corr)),
  x                = hierarchy,
  y                = meanPerf,
  grouping.var = cPairLetters,
  type             = "p",
  bf.message = FALSE
)

## `summarise()` has grouped output by 'SubjectId', 'cPairLetters'. You can
## override using the `.groups` argument.

# Using a mixed logistic model to analyse main effect of Hierarchy 

m1 <- glmer(corr ~ 1 +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2 <- glmer(corr ~ hierarchy +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m1,m2)

	corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.93	1.32 – 2.82	0.001	1.62	1.09 – 2.41	0.017
hierarchy [H2]				1.43	1.15 – 1.77	0.001
Random Effects
σ²	3.29			3.29
τ₀₀	0.59 _SubjectId			0.60 _SubjectId
ICC	0.15			0.15
N	17 _SubjectId			17 _SubjectId
Observations	1632			1632
Marginal R² / Conditional R²	0.000 / 0.151			0.008 / 0.160

anova(m1,m2)

## Data: exp3.inf
## Models:
## m1: corr ~ 1 + (1 | SubjectId)
## m2: corr ~ hierarchy + (1 | SubjectId)
##    npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)   
## m1    2 2010.4 2021.2 -1003.22   2006.4                        
## m2    3 2001.8 2018.0  -997.91   1995.8 10.623  1   0.001117 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# we dont have immediate premise performance but we can lookat the premise X inf interactions

m1 <- glmer(corr ~ 1 +
    (1 | SubjectId), data = exp3_ext, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2 <- glmer(corr ~ pairType +
    (1 | SubjectId), data = exp3_ext, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m3 <- glmer(corr ~ pairType + hierarchy +
    (1 | SubjectId), data = exp3_ext, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m4 <- glmer(corr ~ pairType * hierarchy +
    (1 | SubjectId), data = exp3_ext, family = binomial, control = glmerControl(optimizer = "bobyqa"))
tab_model(m1,m2,m3,m4)

	corr			corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	3.45	2.62 – 4.54	<0.001	1.87	1.38 – 2.53	<0.001	1.63	1.18 – 2.23	0.003	1.58	1.15 – 2.18	0.005
pairType [premise]				4.66	3.86 – 5.61	<0.001	4.68	3.88 – 5.65	<0.001	5.12	3.96 – 6.63	<0.001
hierarchy [H2]							1.32	1.11 – 1.57	0.002	1.41	1.14 – 1.73	0.001
pairType [premise] × hierarchy [H2]										0.82	0.57 – 1.20	0.311
Random Effects
σ²	3.29			3.29			3.29			3.29
τ₀₀	0.30 _SubjectId			0.36 _SubjectId			0.36 _SubjectId			0.36 _SubjectId
ICC	0.08			0.10			0.10			0.10
N	17 _SubjectId			17 _SubjectId			17 _SubjectId			17 _SubjectId
Observations	3264			3264			3264			3264
Marginal R² / Conditional R²	0.000 / 0.084			0.140 / 0.224			0.144 / 0.229			0.143 / 0.228

anova(m1,m2,m3,m4)

## Data: exp3_ext
## Models:
## m1: corr ~ 1 + (1 | SubjectId)
## m2: corr ~ pairType + (1 | SubjectId)
## m3: corr ~ pairType + hierarchy + (1 | SubjectId)
## m4: corr ~ pairType * hierarchy + (1 | SubjectId)
##    npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
## m1    2 3489.1 3501.3 -1742.5   3485.1                           
## m2    3 3195.1 3213.4 -1594.5   3189.1 295.9865  1  < 2.2e-16 ***
## m3    4 3187.0 3211.4 -1589.5   3179.0  10.0826  1   0.001497 ** 
## m4    5 3188.0 3218.4 -1589.0   3178.0   1.0224  1   0.311944    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#create scatterplot
plot(exp3.inf$distance, exp3.inf$jointrank)
#add labels to every point
text(exp3.inf$distance+0.05,exp3.inf$jointrank, labels=exp3.inf$cPairLetters)

# distance
m0.d <- glmer(corr ~ hierarchy  +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m1.d <- glmer(corr ~ hierarchy + distance +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.d <- glmer(corr ~ hierarchy * distance +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m0.d,m1.d,m2.d)

	corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.62	1.09 – 2.41	0.017	0.45	0.26 – 0.80	0.006	0.58	0.30 – 1.14	0.115
hierarchy [H2]	1.43	1.15 – 1.77	0.001	1.44	1.16 – 1.79	0.001	0.85	0.38 – 1.93	0.700
distance				1.62	1.39 – 1.88	<0.001	1.47	1.20 – 1.81	<0.001
hierarchy [H2] × distance							1.22	0.90 – 1.66	0.191
Random Effects
σ²	3.29			3.29			3.29
τ₀₀	0.60 _SubjectId			0.63 _SubjectId			0.63 _SubjectId
ICC	0.15			0.16			0.16
N	17 _SubjectId			17 _SubjectId			17 _SubjectId
Observations	1632			1632			1632
Marginal R² / Conditional R²	0.008 / 0.160			0.040 / 0.193			0.043 / 0.196

anova(m0.d,m1.d,m2.d) #m1.d

## Data: exp3.inf
## Models:
## m0.d: corr ~ hierarchy + (1 | SubjectId)
## m1.d: corr ~ hierarchy + distance + (1 | SubjectId)
## m2.d: corr ~ hierarchy * distance + (1 | SubjectId)
##      npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## m0.d    3 2001.8 2018.0 -997.91   1995.8                          
## m1.d    4 1963.0 1984.5 -977.48   1955.0 40.8532  1  1.641e-10 ***
## m2.d    5 1963.2 1990.2 -976.63   1953.2  1.7079  1     0.1913    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m1.d, type = "pred", terms = c("distance","hierarchy"), axis.title="Inference performance%" )

plot_model(m1.d, type = "diag")

## $SubjectId

## `geom_smooth()` using formula = 'y ~ x'

# jointrank

m0.j <- glmer(corr ~ hierarchy  +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m1.j <- glmer(corr ~ hierarchy + jointrank +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.j <- glmer(corr ~ hierarchy * jointrank +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

tab_model(m0.j,m1.j,m2.j)

	corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.62	1.09 – 2.41	0.017	0.37	0.17 – 0.81	0.012	1.11	0.41 – 3.01	0.844
hierarchy [H2]	1.43	1.15 – 1.77	0.001	1.43	1.16 – 1.78	0.001	0.14	0.04 – 0.56	0.005
jointrank				1.20	1.11 – 1.31	<0.001	1.05	0.94 – 1.18	0.411
hierarchy [H2] × jointrank							1.34	1.13 – 1.58	0.001
Random Effects
σ²	3.29			3.29			3.29
τ₀₀	0.60 _SubjectId			0.61 _SubjectId			0.62 _SubjectId
ICC	0.15			0.16			0.16
N	17 _SubjectId			17 _SubjectId			17 _SubjectId
Observations	1632			1632			1632
Marginal R² / Conditional R²	0.008 / 0.160			0.022 / 0.175			0.033 / 0.186

anova(m0.j,m1.j,m2.j) #m2.j

## Data: exp3.inf
## Models:
## m0.j: corr ~ hierarchy + (1 | SubjectId)
## m1.j: corr ~ hierarchy + jointrank + (1 | SubjectId)
## m2.j: corr ~ hierarchy * jointrank + (1 | SubjectId)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m0.j    3 2001.8 2018.0 -997.91   1995.8                         
## m1.j    4 1984.8 2006.4 -988.42   1976.8 18.972  1  1.327e-05 ***
## m2.j    5 1975.4 2002.4 -982.71   1965.4 11.428  1  0.0007236 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(m2.j, type = "pred", terms = c("jointrank","hierarchy"), axis.title="Inference performance%" )

plot_model(m2.j, type = "diag")

## $SubjectId

## `geom_smooth()` using formula = 'y ~ x'

compare_performance(m1.d,m2.j)

## # Comparison of Model Performance Indices
## 
## Name |    Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## m1.d | glmerMod | 1963.0 (0.998) | 1963.0 (0.998) | 1984.5 (>.999) |      0.193 |      0.040 | 0.160 | 0.447 | 1.000 |    0.581 |      -Inf |           0.002
## m2.j | glmerMod | 1975.4 (0.002) | 1975.5 (0.002) | 2002.4 (<.001) |      0.186 |      0.033 | 0.158 | 0.447 | 1.000 |    0.584 |      -Inf |           0.002

test_performance(m1.d,m2.j)

## Name |    Model |       BF
## --------------------------
## m1.d | glmerMod |         
## m2.j | glmerMod | 1.33e-04
## Each model is compared to m1.d.

plot(compare_performance(m1.d, m2.j))

# jointrank + distance

m0.dj <- glmer(corr ~ hierarchy  +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m1.dj <- glmer(corr ~ hierarchy + jointrank * distance +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))
m2.dj <- glmer(corr ~ hierarchy * jointrank * distance +
    (1 | SubjectId), data = exp3.inf, family = binomial, control = glmerControl(optimizer = "bobyqa"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

tab_model(m0.dj,m1.dj,m2.dj)

	corr			corr			corr
Predictors	Odds Ratios	CI	p	Odds Ratios	CI	p	Odds Ratios	CI	p
(Intercept)	1.62	1.09 – 2.41	0.017	0.15	0.00 – 6.72	0.332	0.30	0.00 – 55.98	0.655
hierarchy [H2]	1.43	1.15 – 1.77	0.001	1.45	1.16 – 1.80	0.001	0.47	0.00 – 932.33	0.844
jointrank				1.15	0.72 – 1.83	0.569	1.08	0.57 – 2.07	0.806
distance				1.38	0.25 – 7.43	0.710	1.67	0.16 – 17.13	0.667
jointrank × distance				1.02	0.83 – 1.26	0.850	0.98	0.74 – 1.32	0.918
hierarchy [H2] × jointrank							1.08	0.42 – 2.79	0.871
hierarchy [H2] × distance							0.59	0.02 – 17.99	0.764
(hierarchy [H2] × jointrank) × distance							1.09	0.71 – 1.68	0.678
Random Effects
σ²	3.29			3.29			3.29
τ₀₀	0.60 _SubjectId			0.64 _SubjectId			0.65 _SubjectId
ICC	0.15			0.16			0.17
N	17 _SubjectId			17 _SubjectId			17 _SubjectId
Observations	1632			1632			1632
Marginal R² / Conditional R²	0.008 / 0.160			0.052 / 0.207			0.064 / 0.218

anova(m0.dj,m1.dj,m2.dj)

## Data: exp3.inf
## Models:
## m0.dj: corr ~ hierarchy + (1 | SubjectId)
## m1.dj: corr ~ hierarchy + jointrank * distance + (1 | SubjectId)
## m2.dj: corr ~ hierarchy * jointrank * distance + (1 | SubjectId)
##       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m0.dj    3 2001.8 2018.0 -997.91   1995.8                         
## m1.dj    6 1948.4 1980.8 -968.21   1936.4 59.386  3  7.952e-13 ***
## m2.dj    9 1941.9 1990.5 -961.97   1923.9 12.489  3   0.005881 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compare_performance(m1.d,m2.j,m2.dj)

## # Comparison of Model Performance Indices
## 
## Name  |    Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## --------------------------------------------------------------------------------------------------------------------------------------------------------------
## m1.d  | glmerMod | 1963.0 (<.001) | 1963.0 (<.001) | 1984.5 (0.952) |      0.193 |      0.040 | 0.160 | 0.447 | 1.000 |    0.581 |      -Inf |           0.002
## m2.j  | glmerMod | 1975.4 (<.001) | 1975.5 (<.001) | 2002.4 (<.001) |      0.186 |      0.033 | 0.158 | 0.447 | 1.000 |    0.584 |      -Inf |           0.002
## m2.dj | glmerMod | 1941.9 (>.999) | 1942.0 (>.999) | 1990.5 (0.048) |      0.218 |      0.064 | 0.165 | 0.442 | 1.000 |    0.571 |      -Inf |           0.002

test_performance(m1.d,m2.j,m2.dj)

## Name  |    Model |       BF
## ---------------------------
## m1.d  | glmerMod |         
## m2.j  | glmerMod | 1.33e-04
## m2.dj | glmerMod |    0.051
## Each model is compared to m1.d.

plot(compare_performance(m1.d,m2.j,m2.dj))

plot_model(m1.dj, type = "pred", terms = c("jointrank","distance"), axis.title="Inference performance%" )

## Could not compute variance-covariance matrix of predictions. No confidence intervals are returned.

plot_model(m1.dj, type = "pred", terms = c("jointrank","distance", "hierarchy"), axis.title="Inference performance%" )

plot_model(m2.dj, type = "pred", terms = c("jointrank","distance", "hierarchy"), axis.title="Inference performance%" )

plot_model(m2.dj, type = "diag")

## $SubjectId

## `geom_smooth()` using formula = 'y ~ x'

Sidetrack: Whats the relationship between jointrank and distance in explaining peoples inference performance regardless of session?

Interim summary

time and sleep-dependent effects are strongly affected by strength of memory encoding
we can see from these 3 datasets that depending on the strength of memory encoding we see time (including overnight sleep) lead to improved or stagnating/diminshed inference performance compared to “baseline”
an understudied aspect of implied serial rank order effects, joint rank of items seems to have a strong effect on inference performance when opportunity to learn is limited
Ellenbogen et al. (2007) showing distance affect without signs of joint rank, similar to Behrens&Bird.. apparently with enough training jointrank effect vanishes which makes sense, with enough training you form structured understanding
- this also highlights the importance of taking into account encoding strength what looking at time/sleep dependent outcomes

	20-min	12-hr	24-hr
B-D	63 (9.5)	74 (7.0)	68 (10.2)
B-E	54 (9.8)	80 (5.6)	86 (6.0)
C-E	40 (11.8)	72 (6.5)	71 (9.4)

	12-hr Wake	12-hr Sleep
B-D	75 (9.6)	73 (10.5)
B-E	69 (8.8)	93 (4.6)
C-E	73 (8.4)	71 (10.7)

Werchan & Gómez (2013) also showing a joint rank of items effect based on the summary statistics (didn’t find sig distance)

	Sleep/reinforcement	Wake/reinforcement	Sleep/no-reinforcement	Wake/no-reinforcement
B > D	76 (7.6)	67 (6.6)	52 (8.6)	46 (10.0)
B > E	85 (5.1)	61 (9.5)	48 (5.9)	46 (6.2)
C > E	60 (8.7)	67 (7.9)	55 (5.7)	58 (8.6

the joint rank of items effect is present even when there is only weak distance effects (e.g: Dataset 1) and interacts with time of encoding consistently (all three datasets)
when strength of memory encoding is low we see that inference performance is driven by low joint rank of items ( Dataset 1, Dataset 2 ), when strength of memory encoding is high ( Dataset 3 ) we see that inference performance is driven by high joint rank of items and highly distant pairs

Learning and Memory

2023-02-13

Version 1: Strength of memory encoding and rank of items affect delayed inference performance

Dataset 1 (Within-subject Transitive Inference with 27h and 3h delayed testing) by Tamas n = 70

Dataset 2 (Between-subject Transitive Inference with 12h delayed testing) by Lorena n=24 (only condition PRE)

No random slope

With random slope

Dataset 3: (Within-subject Transitive Inference with 24h and 0h delayed testing) - Reanalysis of Berens & Bird (2022) n=17

Interim summary

Version 2: Immediate test of inference pairs (Y/N) affects delayed inference performance