Xjyuk

I want to understand how some dependent variable y, depends on a known relationship with independent variable x, but also how x potentially interacts with a high dimensional complex set of features (microbiome data which can be represented as thousands of predictors per observation of y). Therefore, I know that the general form of the model is:

y ~ m1*x + b

However I would like to add an interaction term, generating the model:

y ~ m1*x + m2*x*microbiome + b

Where microbiome is the very large microbiome feature set.

I think the microbiome data essentially interacts with the independent variable x to affect y, but I don't know how. There are thousands of species, and my guess is that certain combinations will be predictive of this interaction and explain a lot of variation, but a priori I don't know which. I also suspect that the microbiome features relate to each other in a non-linear way. Essentially I want to use a machine learning approach to figure that out. I am aware of "boosted" regressions, where you use a machine learning algorithm on the residuals, however I want to specify something a bit more mechanistic than that.

If anyone could suggest a method to do something like this (especially if it can be implemented in R) I would be very interested. If there is a way to use the residuals of the model to do this, I would also be interested. It's worth noting that in the actual application I have many more predictors in the model, many of which have non-linear relationships with y.

The easiest way to implement model like

y ~ m1*x + m2*x*microbiome + b

would be to replace microbiome with a dense neural network

y ~ m1*x + m2*x*nn(microbiome) + b

so that the neural network nn would reduce dimensionality (to single or multiple dimensions, depending of number of units in the output layer) and do the feature engineering for you. The nice part is that it would let you to keep the assumed form of the model, but the neural network would deal with the extra features for you.

This can be easily done in frameworks like Keras, that are designed to deal with large datasets and scale nicely. In Keras, this would translate to something like the model definition below. To understand the code, you would probably need to dive deeper into Keras, but hopefully many tutorials are available online.

from keras.models import Model
from keras.layers import Input, Dense, multiply, concatenate

x_inp = Input(shape=(1,))
microbiome_inp = Input(shape=(k,))

# 3-layer neural network
nn = Dense(200, activation='relu')(microbiome_inp)
nn = Dense(50, activation='relu')(nn)
nn = Dense(1)(nn)

# x*nn(microbiome)
mul = multiply([x_inp, nn])

# m1*x + m2*x*nn(microbiome) + b
conc = concatenate([x_inp, mul])
out = Dense(1)(conc)

model = Model(inputs=[x_inp, microbiome_inp], outputs=out)
model.compile(optimizer='adam', loss='mean_squared_error')
model.fit([x, microbiome], y)