modal-visualization/LinearModels_func.py at main · PhilLaboratory/modal-visualization · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
# Functions for fitting linear models and computing partial R-squares

import pandas as pd
import statsmodels.formula.api as smf
from scipy.stats import chi2

#%%
# Goodness/Poss X Time interaction
# Fit linear model, do model comparison with reduced models, and compute partial r-squares for interaction terms
def lm_compare(modal, df):
  full_interaction = True
  random_intercepts = False
  df_select = df[['trialNo', 'Factor_1', 'Factor_2', f'{modal}_resp_sp', f'{modal}_resp_rf']]

  df_long = pd.melt(df_select, id_vars=['trialNo', 'Factor_1', 'Factor_2'], value_vars=[f'{modal}_resp_sp', f'{modal}_resp_rf'],
                  var_name='time_cond', value_name='response')

  # model formulas
  if full_interaction:
    f0 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_1:time_cond + Factor_2:time_cond + Factor_1:Factor_2:time_cond + C(trialNo)'
    f1 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_2:time_cond + Factor_1:Factor_2:time_cond + C(trialNo)' # no f1*time
    f2 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_1:time_cond + Factor_1:Factor_2:time_cond + C(trialNo)' # no f2*time
  elif not full_interaction:
    f0 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_1:Factor_2 + Factor_1:time_cond + Factor_2:time_cond'
    f1 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_1:Factor_2 + Factor_2:time_cond' # no f1*time
    f2 = 'response ~ Factor_1 + Factor_2 + time_cond + Factor_1:Factor_2 + Factor_1:time_cond' # no f2*time

  # specify models
  if not random_intercepts:
    # linear model without random intercepts
    model_full = smf.ols(f0, data=df_long).fit()
    model_r1 = smf.ols(f1, data=df_long).fit()
    model_r2 = smf.ols(f2, data=df_long).fit()
  elif random_intercepts:
    model_full = smf.mixedlm(f0, data=df_long, groups=df_long['trialNo']).fit(method="nm")
    model_r1 = smf.mixedlm(f1, data=df_long, groups=df_long['trialNo']).fit(method="nm")
    model_r2 = smf.mixedlm(f2, data=df_long, groups=df_long['trialNo']).fit(method="nm")

  # likelihood ratio test
  # full model vs. reduced model 1
  lr_stat1, p_value1 = get_llr_stats(model_full, model_r1)
  # full model vs. raduced model 2
  lr_stat2, p_value2 = get_llr_stats(model_full, model_r2)

  # compute partial r-sq for interaction terms
  partial_r2_f1xtime = get_partial_r2(model_full, model_r1)
  partial_r2_f2xtime = get_partial_r2(model_full, model_r2)

  # note the direction of the change in effect
  f1_greater = 'sp'
  f2_greater = 'sp'
  if model_full.params[f'Factor_1:time_cond[T.{modal}_resp_sp]'] <= 0: #<- effect of factor lower in sp than rf
    f1_greater = 'rf'
  if model_full.params[f'Factor_2:time_cond[T.{modal}_resp_sp]'] <= 0:
    f2_greater = 'rf'

  result = {"modal judgment": modal,
            "chi-sq1": lr_stat1,
            "p_val1": p_value1,
            "chi-sq2": lr_stat2,
            "p_val2": p_value2,
            "full model r-sq": model_full.rsquared,
            "partial_r-sq_f1xtime": partial_r2_f1xtime,
            "partial_r-sq_f2xtime": partial_r2_f2xtime,
            "f1 effect is greater in": f1_greater,
            "f2 effect is greater in": f2_greater}

  return result

#%%
# Goodness/Poss X Modal term (pairwise)
# Fit linear model, do model comparison with reduced models, and compute partial r-squares for interaction terms
def factor_modal_interation(df, modal1, modal2, time='rf'):

  df_select = df[['trialNo', 'Factor_1', 'Factor_2', f'{modal1}_resp_{time}', f'{modal2}_resp_{time}']]
  df_select.columns = ['trialNo', 'Factor_1', 'Factor_2', f'0_{modal1}_resp_{time}', f'1_{modal2}_resp_{time}']

  df_long = pd.melt(df_select, id_vars=['trialNo', 'Factor_1', 'Factor_2'], value_vars=[f'0_{modal1}_resp_{time}', f'1_{modal2}_resp_{time}'],
                  var_name='modal_term', value_name='response')

  # model formulas
  f0 = 'response ~ Factor_1 + Factor_2 + modal_term + Factor_1:modal_term + Factor_2:modal_term + Factor_1:Factor_2:modal_term + C(trialNo)'
  f1 = 'response ~ Factor_1 + Factor_2 + modal_term + Factor_2:modal_term + Factor_1:Factor_2:modal_term + C(trialNo)' # no f1*time
  f2 = 'response ~ Factor_1 + Factor_2 + modal_term + Factor_1:modal_term + Factor_1:Factor_2:modal_term  + C(trialNo)' # no f2*time

  # specify models
  model_full = smf.ols(f0, data=df_long).fit()
  model_r1 = smf.ols(f1, data=df_long).fit()
  model_r2 = smf.ols(f2, data=df_long).fit()

  # compute partial r-square
  partial_r2_f1xmodal = get_partial_r2(model_full, model_r1)
  partial_r2_f2xmodal = get_partial_r2(model_full, model_r2)

  # model comparison/likelihood ratio test
  lr_stat1, p_value1 = get_llr_stats(model_full, model_r1)
  lr_stat2, p_value2 = get_llr_stats(model_full, model_r2)

  # note the direction of the change in effect
  f1_greater = modal2
  f2_greater = modal2
  if model_full.params[f'Factor_1:modal_term[T.1_{modal2}_resp_rf]'] <= 0: #<- effect of factor lower in modal2
    f1_greater = modal1
  if model_full.params[f'Factor_2:modal_term[T.1_{modal2}_resp_rf]'] <= 0:
    f2_greater = modal1

  result = {"modal1": modal1,
            "modal2": modal2,
            "chi-sq1": lr_stat1,
            "p_val1": p_value1,
            "chi-sq2": lr_stat2,
            "p_val2": p_value2,
            "full model r-sq": model_full.rsquared,
            "partial_r-sq_f1xmodal": partial_r2_f1xmodal,
            "partial_r-sq_f2xmodal": partial_r2_f2xmodal,
            "f1 effect is greater in": f1_greater,
            "f2 effect is greater in": f2_greater}

  return result

#%%

# Compute partial r-square for uinque term in the full model
def get_partial_r2(model_full, model_reduced):
  # Full model R-square
  r2_full = model_full.rsquared
  # Reduced model R-square
  r2_reduced = model_reduced.rsquared
  # Partial R-square for the unique term
  partial_r2 = (r2_full - r2_reduced) / (1 - r2_reduced)

  return partial_r2

# Model comparison using likelihood ratio test, return p-value
def get_llr_stats(model_full, model_reduced):
  lr_stat = -2 * (model_reduced.llf - model_full.llf)
  df_diff = model_full.df_model - model_reduced.df_model
  p_value = chi2.sf(lr_stat, df_diff)
  return lr_stat, p_value