The following question is sort of a concrete adaption of a post on StatsExchange [https://stats.stackexchange.com/questions/10182/intraclass-correlation-coefficient-vs-f-test-one-way-anova/11732#11732].

The following R script runs just fine:

library(reshape)

J1 <- c(9,6,8,7,10,6)
J2 <- c(2,1,4,1,5,2)
J3 <- c(5,3,6,2,6,4)
J4 <- c(8,2,8,6,9,7)
Subject <- c('S1', 'S2', 'S3', 'S4', 'S5', 'S6')
sf <- data.frame(Subject, J1, J2, J3, J4)
sf.df <- melt(sf, varnames=c("Subject", "Rater"))
anova(lm(value ~ Subject, sf.df))
anova(lm(value ~ Subject*variable, sf.df))

However, when I try to translate this into Python, I run into an exception that gets thrown from scipy.

Here’s the equivalent Python script:

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# Create the data
J1 = [9, 6, 8, 7, 10, 6]
J2 = [2, 1, 4, 1, 5, 2]
J3 = [5, 3, 6, 2, 6, 4]
J4 = [8, 2, 8, 6, 9, 7]
Subject = ['S1', 'S2', 'S3', 'S4', 'S5', 'S6']

# Create the DataFrame
sf = pd.DataFrame({'Subject': Subject, 'J1': J1, 'J2': J2, 'J3': J3, 'J4': J4})
print(sf)

# Melt the DataFrame
sf_melted = pd.melt(sf, id_vars=['Subject'], var_name='Rater', value_name='Value')
print(sf_melted)

# Perform ANOVA
model1 = ols('Value ~ Subject', data=sf_melted).fit()
anova_table1 = sm.stats.anova_lm(model1, typ=2)
print(anova_table1)

model2 = ols('Value ~ Subject * Rater', data=sf_melted).fit()
anova_table2 = sm.stats.anova_lm(model2, typ=2)
print(anova_table2)

The ANOVA calculation for the first model runs just fine, but there seems to be an issue with the second model.

Here’s the full stack trace:

/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/statsmodels/regression/linear_model.py:1717: RuntimeWarning: divide by zero encountered in double_scalars
  return np.dot(wresid, wresid) / self.df_resid

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[246], line 26
     23 print(anova_table1)
     25 model2 = ols('Value ~ Subject * Rater', data=sf_melted).fit()
---> 26 anova_table2 = sm.stats.anova_lm(model2, typ=2)
     27 print(anova_table2)

File ~/anaconda3/envs/python3/lib/python3.10/site-packages/statsmodels/stats/anova.py:353, in anova_lm(*args, **kwargs)
    351 if len(args) == 1:
    352     model = args[0]
--> 353     return anova_single(model, **kwargs)
    355 if typ not in [1, "I"]:
    356     raise ValueError("Multiple models only supported for type I. "
    357                      "Got type %s" % str(typ))

File ~/anaconda3/envs/python3/lib/python3.10/site-packages/statsmodels/stats/anova.py:84, in anova_single(model, **kwargs)
     81     return anova1_lm_single(model, endog, exog, nobs, design_info, table,
     82                             n_rows, test, pr_test, robust)
     83 elif typ in [2, "II"]:
---> 84     return anova2_lm_single(model, design_info, n_rows, test, pr_test,
     85                             robust)
     86 elif typ in [3, "III"]:
     87     return anova3_lm_single(model, design_info, n_rows, test, pr_test,
     88                             robust)

File ~/anaconda3/envs/python3/lib/python3.10/site-packages/statsmodels/stats/anova.py:207, in anova2_lm_single(model, design_info, n_rows, test, pr_test, robust)
    205 LVL = np.dot(np.dot(L1,robust_cov),L2.T)
    206 from scipy import linalg
--> 207 orth_compl,_ = linalg.qr(LVL)
    208 r = L1.shape[0] - L2.shape[0]
    209 # L1|2
    210 # use the non-unique orthogonal completion since L12 is rank r

File ~/anaconda3/envs/python3/lib/python3.10/site-packages/scipy/linalg/_decomp_qr.py:129, in qr(a, overwrite_a, lwork, mode, pivoting, check_finite)
    125     raise ValueError("Mode argument should be one of ['full', 'r',"
    126                      "'economic', 'raw']")
    128 if check_finite:
--> 129     a1 = numpy.asarray_chkfinite(a)
    130 else:
    131     a1 = numpy.asarray(a)

File ~/anaconda3/envs/python3/lib/python3.10/site-packages/numpy/lib/function_base.py:603, in asarray_chkfinite(a, dtype, order)
    601 a = asarray(a, dtype=dtype, order=order)
    602 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
--> 603     raise ValueError(
    604         "array must not contain infs or NaNs")
    605 return a

ValueError: array must not contain infs or NaNs

Things I’ve tried:

1. Changing library versions. I was originally on numpy 1.22.4 and scipy 1.12, and now I’m on numpy 2.0.0 and scipy 1.13.1.
2. I’ve tried enforcing categorical variables in the linear model by writing Value ~ Subject * C(Rater) for the second model. I wonder if it’s not treating the data as categorical….but from what I’ve read statsmodels is supposed to correctly infer that any string variable is categorical, so this may be the wrong idea.

Can any R/Python/statsmodels experts help me resolve this?