How to perform a regression with categorical variables as input and numerical output? I am a chemist and I am currently working on a machine learning/statistical learning based project. This question is related to a previous question that I asked before, but is more about the learning method.
I have a set of molecules, and I have to predict the spectral peak for each molecule. That's the chemistry part. The statistical modelling part is that I have an array for each molecule as input, which have binary data that looks like this:
molecule1_fingerprint = [0,1,0,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,...,0,0,0,0,1]

Each number indicates the presence or absence of a certain feature in the molecule. I can also use another type of numbers (called hashed fingerprints):
molecule1_hashed = [0,2,0,6,1,0,7,15,...,0,7,0]

Here similar features are expressed by one number so some adjacent numbers in the original array gets added together, reducing the size of the input.
So my input data would look something like this:
X_train = [[0,0,0,1,0,0,1,...,0,0,1,0],
           [0,1,0,0,1,0,1,...,0,0,1,0],
           : : : : : : : : : : : : : :
           : : : : : : : : : : : : : :
           [0,1,0,0,1,1,1,...,0,1,0,0]]

My target is would be a floating point number, (I could use integers if required):
y_train = [100.6, 20.8, -12.9, 10.8, -19.2, 5.7, ..., 8.9, -40.3]

I am not sure how to approach this problem. I have tried using support vector regression, but it gives very poor results (score of -0.003). My intuition from chemistry is that the target y value would have contributions from each of the feature in the input, (i.e. the presence of certain feature would add or subtract a certain amount, and the total value would be obtained by adding those contributions.)
Another problem is that I cannot standardize the y values, because I need their exact values to be predicted. So for example, if I standardize for one dataset, there is no guarantee that the standardization would be the same for another dataset. I am using scikit-learn for my calculation.
I have uploaded a sample dataset that I am working with, on this github page. The arrays can be loaded with numpy.loadtxt(), and are in the scikit-learn orientation. There are 65 molecules in the dataset, and 732 features. The first 731 features are hashed fingerprints, the last feature can be ignored for now.
Any advice about how to choose a method and how to handle the data?

Your data consist of nearly 10x more variables than observations, with 65 rows and nearly 732 columns.  This is going to be tough, but we will persevere.
Because the data have so many variables, you're going to need a method which can deal with a lot of variables and avoid overfitting.  There are a few ways to do this, but first we need a baseline measure of performance.  Sklearn's DummyRegressor is a good place to start.
import numpy as np

from sklearn.dummy import DummyRegressor
from sklearn.metrics import mean_squared_error, make_scorer
from sklearn.model_selection import GridSearchCV, cross_validate
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestRegressor

X = np.loadtxt('https://raw.githubusercontent.com/shoubhikraj/molecular-modelling/main/uv-vis-pred/X_train.txt')
y = np.loadtxt('https://raw.githubusercontent.com/shoubhikraj/molecular-modelling/main/uv-vis-pred/y_train.txt')

rmse = make_scorer(mean_squared_error, squared=False)
results = cross_validate(DummyRegressor(), X, y, scoring = rmse)
results['test_score'].mean()
>>>46.09

Any model with an RMSE above 46 is not worth consideration, because it seems we can a achieve a lower RMSE just by guessing the sample mean.  Let's move on to a linear model.
Because we have so many variables, the linear model must either a) project the variables onto a lower dimensional space, and or b) use regularization.  Let's use both in the following Lasso model
model = Pipeline([
    ('scale', StandardScaler()),
    ('pca',PCA()),
    ('lm', Lasso())
])

param_grid = {
    'pca__n_components': [3, 5, 7, 9, 11, 40],
    'lm__alpha': np.logspace(-3, 3, 5)
}

gscv = GridSearchCV(model, param_grid=param_grid, cv=10, scoring = rmse)

results = cross_validate(gscv, X, y, cv=10, scoring = rmse, verbose=4)

results['test_score'].mean()
>>>35.75

Awesome, through nested cross validation, the model I've used is estimated to have an out of sample RMSE of 35.75.  Its worth considering another model, the random forest. This model is great when we have lots of variables because it randomly selects a subset to create splits.  Here are the cross validation results.
model = RandomForestRegressor(n_estimators=5000, max_features=1)
results = cross_validate(model, X, y, cv=10,  scoring = rmse)
results['test_score'].mean()
>>>37

Promising results.  I bet with a grid search cross validation we could do better.  Additionally, I would perhaps try boosted trees as they seem to do very well.
Hope this helps.
