PyMC3 posterior prediction example; how to incorporate a matrix of betas?

Thanks! That did it. The problem was:

# specify model
with pm.Model(coords=coords) as sk_lm:
    # data
    feature_data = pm.Data('feature_data', features, dims=('obs_id', 'features'))

    # priors
    alpha = pm.Normal('alpha', mu=0, sigma=1)
    betas = pm.Normal('betas', mu=[40, 90, 50], sigma=5, dims='features')

    # model error
    sigma = pm.Exponential("sigma", lam=1)

    # matrix-dot products
    m1 = pm.math.matrix_dot(feature_data, betas)  # here

I was using the literal features data instead of the pm.Data object.