OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Translated R Script Throws Error in Python

  • Thread starter Thread starter hobscrk777
  • Start date Start date
H

hobscrk777

Guest
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:

Code:
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:

Code:
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:

Code:
/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?
<p>The following question is sort of a concrete adaption of a post on StatsExchange<a href="https://stats.stackexchange.com/que...efficient-vs-f-test-one-way-anova/11732#11732"> [https://stats.stackexchange.com/que...efficient-vs-f-test-one-way-anova/11732#11732]</a>.</p>
<p>The following R script runs just fine:</p>
<pre><code>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))
</code></pre>
<p>However, when I try to translate this into Python, I run into an exception that gets thrown from scipy.</p>
<p>Here's the equivalent Python script:</p>
<pre><code>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)
</code></pre>
<p>The ANOVA calculation for the <em>first</em> model runs just fine, but there seems to be an issue with the second model.</p>
<p>Here's the full stack trace:</p>
<pre><code>/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
</code></pre>
<p>Things I've tried:</p>
<ol>
<li>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.</li>
<li>I've tried enforcing categorical variables in the linear model by writing <code>Value ~ Subject * C(Rater)</code> 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.</li>
</ol>
<p>Can any R/Python/statsmodels experts help me resolve this?</p>
 

Latest posts

Top