Statsmodels：使用无法识别的数据结构编写公式

Question

我正在尝试进行多元回归分析。 这里显示我的脚本：

import numpy  as np
import xarray as xr
import cftime
import statsmodels.api as sm

# load two 1D SSTA data
nc_file = "nino3_hadley_197901_202212_monthly.nc" 
ds = xr.open_dataset(nc_file)
n3 = ds.n34

nc_file = "nino4_hadley_197901_202212_monthly.nc" 
ds = xr.open_dataset(nc_file)
n4 = ds.n34

# Load the 3D precipitation data
nc_file2 = "precip.mon.mean.197901-202212.nc"
ds = xr.open_dataset(nc_file2)
pr = ds.precip

# Calculate the climatology
pr_climatology = pr.groupby('time.month').mean('time')

# Calculate the anomaly
pr_anomaly = pr.groupby('time.month') - pr_climatology

# multiple regression
X = np.column_stack((n3, n4))
# add the constant value of 1 to the leftside array. 
X = sm.add_constant(X)

一切正常，直到我运行如下所示的命令：

# Perform the regression
model = sm.OLS(pr_anomaly[:,36,72], X).fit()

错误显示： ValueError：无法识别的数据结构： /

图片：

有人可以帮忙吗？我是 python 的新手。

print(n3)
<xarray.DataArray 'n34' (time: 528)>
array([-0.04474 , -0.164384, -0.060881, ..., -1.03504 , -0.729144, -0.845494],
      dtype=float32)

print(pr)
<xarray.DataArray 'precip' (time: 528, lat: 72, lon: 144)>
[5474304 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 1979-01-01 1979-02-01 ... 2022-12-01
  * lat      (lat) float32 -88.75 -86.25 -83.75 -81.25 ... 83.75 86.25 88.75
  * lon      (lon) float32 1.25 3.75 6.25 8.75 11.25 ... 351.2 353.8 356.2 358.8