Skip to content

My words, sometimes technical, sometimes not

Taking notes is life changing.

I constantly feel that my memory is failing me more and more even for things that I put effort of retaining.
Why is a more complex topic, I'm not sure, some medical condition? chronical bad sleep? nothing wrong just this is it and everyone has this poor memory?

The feeling is like when you read some code you wrote few weeks/months ago and you can't remember anything nor why you did it this way. Well, that but for a lot of stuff.

Despite the cause I've decided to help myself and stop thinking "I will remember this, it's super important" and I began taking notes regularly. This helped me in a huge way for personal topics but also for work and technical subjects.

It helps in two ways:

1) Writing helps memorizing things, you are going through it again and you need to organize a bit your ideas to put them on paper.
2) I have a place to rest on where I can find information I should have memorized.

I started doing it randomly for everyday things, like things I want to buy for the house, places I would like to visit, places I've gone.
This Christmas I'm also writing the gifts I gave, and what I received. I hate the feeling of someone telling me "Hey, this is the shirt you gave me!" and I'm speechless for a second and it's obvious I don't remember, or to avoid giving twice the same gift.
I'm expanding these notes as I encounter situations where I regret the status of my memory.

But taking notes is also huge for me in terms of Work and research.

Some examples:

  • Daily notes on what I've done in my full time job (YYYY-MM-DD). I can always come back if someone asks, sometimes I put some briefs about decisions made on meetings. Specially useful for Mondays where I forget what I did/talked on Friday.
  • Notes on how I did something that I did for the first time. In my job or in my side businesses, including decisions, commands, installation processes. I find myself coming back to these really often and I think this saves my huge amounts of time, stress and frustration. It's surprising how I forget THE command I've been running twice a day for some deployment test if I have to come back to that topic two weeks later. Something that seems impossible to forget at some point just vanishes. Also, super useful for my side business where I setup something once (DB, platform, server, whatever) and I need to revisit how I did it or what resource I followed in case I have to do it again.
  • Studying. I'm also taking side notes on books and papers that I later write down in the form of more structured texts. What I understood about papers, doubts, summaries. As usual, while reading I feel I understand but then if I need to revisit the paper without notes it feels like I need to re read it. Many times I read my notes and I feel surprised as I don't remembered a lot of contents, is like the note was given to me by someone else.

Motivation

It's a tough topic because everyone has a different perception and it's hard to share emotions but I feel what I call motivation really has an impact on me. Sometimes is a silent impact, I'm talking now about lack ok motivation, noticeable only when it comes back.
You don't know what you don't know. You don't know how much energy and willingness to do things you should have if you never realized it.

My biggest adulthood lack of motivation came from jobs. I can't complain, I work in software/Machine learning, good salaries, remote work if desired, etc, but still, a day to day that is not fulfilling has devastating effects. Even more, probably fully remote work post pandemic plays a role in this despite I think the net effect is positive as I can't value enough the commute time regained.

The difference between being motivated to do something and being indifferent is the largest gap possible. If you truly desire something and work for it the odds of succeding increase exponentially.

Things I've realized that affect my motivation:

  • Sleep. I have regular issues sleeping. Usually after 2 days of bad sleep my brain doesn't work anymore. I work poorly, I get anxious because I can't get anything done nor I can think properly. I sleep bad again and so on. Disaster.
  • Anxiety. It's related to sleep but not always. Feeling I can't do something or that I don't understand some topic or how to solve a problem, makes me anxious and after a while I want to quit. I trust I can understand or solve complex machine learning topics but usually not having a clear North is really stressful. If only I could have someone guiding a bit, everything looks better. Not knowing if you are on the right path and you have a short time to figure out something is horrible. For instance, we need to solve X problem. You ask a few days for research, you accumulate papers or blogposts, but sometimes none is exactly your solution, and then all have something different, and then you don't know what to read and all are long pieces of text, and they mention other topics you don't know. I feel anxious just by writing this, but one needs to make peace with the fact that you can't have a phD in every topic to solve. Maybe in none if you are just in industry and you are responsible for multiple projects, and things will be done suboptimally probably. But it's hard to really accept it and on top of that, is you have a boss or some coleague asking for this, the stress is twice as big, at least for me. Fear of disappointing. Maybe that's another whole note.
  • Daily meetings. Daily standup or just status meetings are terrible for me. I feel the pressure to explain myself every day, to tell someone else what did I do for 8 hours and sometimes, when motivations it's at a low level, is not that much, or maybe I am just researching without having grasped full understanding so it seems unfruitful, and this goes back to anxiety. Every day feeling the need to have something meaningful to share. Currently I don't have daily meetings, I have some freedom to work and then put up to speed my manager. I feel more free and funny enough, I feel that I do more progress on a daily basis this way.
  • Social media. You can see unrealistic success stories all day long that feel like someone found the way to live in a 5 minute revelation and minimum effort. All fake, all lies but you can quickly fall in that trap of feeling less than others. You just have to compare yourself with previous you, not others that are even probably fake. On top of that, when you realize all the time you get back you feel like your days are longer and you can fight the anxiety better. At some point I realized this and when I automatically reached my phone, which happened really often, to open Twitter, I would hit my wrist (without hurting!) and took a deep breath, then drop the phone. You get back 5/10 minutes of time you would use in social media, and that accumulates fast.

I've quit because of this, because of not being happy, because of being stressed every morning by lack of progress. Maybe I wasn't up to the task, maybe if you know exactly what you are doing you can move forward steadily without worrynig about this but I'm sure there was another time in a project where I felt all this lack of motivation and it was not because of "skill issue" but because of continuous urgency. Every morning I could (or maybe not, you don't know until you check your phone) with some urgent message of modifying something because yadda yadda the client. Sometimes it made sense, others.. not so much, and that feeling of being alert and thinking about slack messages all the time was awful, you can't focus, you can't feel you are working in peace in something.

Consequence

Eventually I quit all that and started learning about how to set up a job board, reading , drinking coffee, I took a 6 month sabbatical, at least, industry sabbatical, but I would spend time in the computer learning, working in the job board, etc. But I was free and owner of my own time, no one expecting anything. If I wanted to do nothing and just play with my dog all afternoon, there was no daily standup the next morning to dread. I could enjoy things without worrying.

Unfortunately in those six months I didn't make a fortune nor won the jackpot so I went back to industry but without ruhsing. I looked for jobs that at least, before joining, made me want to join. Maybe because of the area, maybe because of the technology and in some case, maybe because it seemed chill. I had my share of rejections and processes that went nowhere but I ended up in a startup with a really enjoyable day to day structure. I'm really grateful about it, I have a couple of meetings per weeks and that's it. I can do things at my own pace (given some realistic timelines), I can share ideas with my manager but I don't need to impress anyone on a daily basis. You could think that I'm chilling but not, I work hard, I learn a lot of stuff and I put the best of me because I feel comfortable and grateful to how my manager handles the startup.

Conclusions

  • Fix your sleep. Try everything and don't give up for too much time. You can try:
    • Magnesium
    • Having dinner earlier
    • Start relaxing one hour before (lights down, no screens, reading a bit, etc). This is the most unrealistic for me, can be done once, not regularly if you are anxious.
    • Drink some relaxing tea.
    • Legal drugs (alprazolam or similar, check with doctor.)
    • Nose strips
    • Meditate

You need to find what's the main cause but if it's only "being anxious", start trying stuff, at least you will feel doing something about it.

  • If you can afford it, change your job, take a sabbatical, take some time and look for what you really want.

Expectation Maximization for Poisson process

The problem

What if you want to run a regression to estimate the coefficients and retrieve the data generating process of some data BUT you are in a scenario with censored data? How does that affect my estimation? Can I do better?

Censored data

In short, you have censored data when some observation is unobserved or constrained because of some specific reason which is natural and unavoidable.

Unobserved

For instance, in survival analysis, if the subject is not yet dead (luckily) you can not observe the time of death yet and hence you don't know if he will die tomorrow or in two years. You can't observe yet the death date.

Constrained

For goods with limited demand, once you deplete your stock, by definition you can't sell more and you can't observe all the demand that you would have if more items were able to be sold. Or an Airline , when selling seats, can't observe the full demand after all the seats have been sold, maybe more people were willing to flight.

We will simulate an example of the latest case, for an airline with simple toy data.

import numpy as np
import statsmodels.api as sm
from scipy.stats import poisson
import matplotlib.pyplot as plt

We create some fake data, where seats sold are generated from a poisson distribution, with an intercept (2) and a coefficient related to how many days are left to the departure (-0.2), the further from the departure, the less seats are sold.

What we also do, is that on the last date, we constrain the demand. We don't get the full seats that would have been sold under a poisson distribution but the amount minus 3. Only ourselves know that because we want to tweak the data. In real life we would see data similar to the one generated and we couldn't see the actual sales that would have been present in a complete poisson generating process (that actually is how the tickets are sold in our example)

# Step 1: Generate data
np.random.seed(100)
n_samples = 20
n_features = 1
substract_sales = 3

# Generate features
X = np.random.normal(size=(n_samples, n_features))
X = np.array(range(n_samples+1, 1, -1))
X = sm.add_constant(X)  # Add intercept term

# True coefficients
beta_true = np.array([2, -0.2])

# Generate Poisson-distributed target values
linear_pred = np.dot(X, beta_true)
lambda_true = np.exp(linear_pred)
y = np.random.poisson(lambda_true)
# y

# Right-censoring threshold in the last date (when the fare is closed)
censor_threshold = np.concatenate([(y+1)[:-1],np.array([y[-1]-substract_sales])])
y_censored = np.where(y > censor_threshold, censor_threshold, y)
is_censored = (y > censor_threshold)
Maybe it's easier to visualize.
The green curve is the actual lambda parameter for each data, based on the intercept and the real coefficient multiplied by days to departure.
The blue dots are the actual data points we see for all dates except from the last one.
The grey dot is the data we should see if there was no limit of how many seats to sell, following the poisson distribution.
The red dot is the actual data point we have, the last day the airline sold all the remaining seats and couldn't fulfill the true demand and hence we see lower sales than the ones generated by the poisson distribution. 

# Plot the data
plt.xlabel('Days to departure')
plt.ylabel('Bookings')
plt.title('Data Plot')
plt.plot(X[:-1, 1], y[:-1], 'o')
plt.plot(X[:, 1], lambda_true, label='Lambda True', color="green")
# plt.plot(X[:, 1], lambda_est, label='Lambda Estimated')
plt.scatter(X[-1, 1], y[-1:],  color="grey", label="Real unseen demand" )
plt.scatter(X[-1, 1], y[-1:]-substract_sales, color='red', label='Constrained')
plt.legend()
plt.show()

Image

Esimation
Poisson Regression

First we estimate the parameters using a poisson regression. We assume the data follows a poisson process and we don't do anything to manage the constrained data.
The results are kind of ok, we get some closeness to the true parameters but it is clearly biased.

model_censored = sm.GLM(y_censored, X, family=sm.families.Poisson())
results_censored = model_censored.fit()
print(f'Estimated Censored coefficients: {results_censored.params}')

## Estimated Censored coefficients: [ 1.86976922 -0.21124541]
Expectation Maximization

OK, here is our alternative. What we can do, a bit more sophisticated, is to apply EM algorithm, which is a iterative process, to estimate the coefficients of the poisson regression, including the knowledge we have, that there might be constrained data.

# Step 2: Initialize parameters
beta_est = np.zeros(n_features + 1)
tol = 1e-4
max_iter = 1000

for iteration in range(max_iter):
    # Step 3: E-step
    # Estimate the expected values of censored data
    lambda_est = np.exp(np.dot(X, beta_est))
    expected_y = np.where(
        is_censored,
        (censor_threshold + 1) / (1 - poisson.cdf(censor_threshold, lambda_est)),
        y_censored
    )
    # Ensure expected_y values are valid
    expected_y = np.nan_to_num(expected_y, nan=np.mean(y_censored), posinf=np.max(y_censored), neginf=0)
    # Step 4: M-step
    # Update parameter estimates using Poisson regression
    model = sm.GLM(expected_y, X, family=sm.families.Poisson())
    results = model.fit(method="lbfgs")
    # results = model.fit_regularized(L1_wt=0, alpha=0.1) 
    new_beta_est = results.params

    # Check convergence
    if np.linalg.norm(new_beta_est - beta_est) < tol:
        break

    beta_est = new_beta_est

## <string>:7: RuntimeWarning: divide by zero encountered in divide

print(f'Estimated coefficients: {beta_est}')


## Estimated coefficients: [ 2.11015996 -0.23573144]
We can see our estimates are of course not perfect but the intercept is closer to the true parameter.
Let's visualize the results to have a clear intuition.

uncensored_predictions = results.predict(X)
censored_predictions = results_censored.predict(X)
Promising!
For the furthest dates we see a bias in both approaches, we are not doing better than the Poisson regression, but neither worse.
As we move closer the the departure, the EM algorithm predictions get closer and closer to the true lambdas, while the regular poisson regression continues it's biased trajectory.
The EM procedure adjusted the coefficients to better match the constrained data.

plt.plot(X[:, 1], lambda_true, label='Lambda True', color="green")
plt.plot(X[:, 1], uncensored_predictions, label='Uncensored Predictions (EM)', color="blue")
plt.plot(X[:, 1], censored_predictions, label='Censored Predictions', color="red")
plt.xlabel('Days to departure')
plt.ylabel('Lambda - Poisson regression')
plt.title('Lambda Predictions')
plt.legend()
plt.show()

Image

Reflections on quitting my ML job

As I'm starting my sabbatical journey I am reading some posts of Jason Liu since he seems to have a career similar to something I would like to if I got into ML and LLM riding the hype wave. Furthermore I find his writing pleasant and his content could be useful even as a solopreneur.

One of my goals for this period of time is to write more. To take more notes since I struggle memorizing without them (or even with them but at least I can read a summary later) and I want to actually do more. More of everything, less thinking about and more actual doing. From shipping products, to really learning and that includes writing.

Reading his post and starting now my journey, a real action would be to take notes about his post, which I found useful, at least some bits of it. Not taking notes would be less writing, less doing, less memorizing.

Choosing

Right now, I feel more at this point

"This despair arises from the realization of one's absolute freedom and the responsibility for creating one's own essence and purpose."

Lately, last few years probably, I have started to realize that quote. We can go a long way by following the usual paths, at least usual to our surroundings. In my case, without major distress that turns everything upside down, it was high school, university, get a job, get married, retire. At least, that's something I (many?) thought as given, a fate that if didn't screw it, it would happen automatically.

And as I was applied and I did good in school, I got the first 3 steps quite easily. I never really thought of going out of that path and saw the ones doing so as outliers or with greater safety nets to fail. In some cases it could be true, in others it was probably me being short sighted, I guess it was not my fault but just lack of adult life.

As years go by, I start questioning the meaning of what I'm doing. Is my job meaningful? Does it create value or help anyone? Corporate world feels like a charade. Despite I was working in analytics and studying and building ML models, I couldn't see if that was worth the effort and time. Once you go past the hype of learning models and cool tech bits, looking over that it feels empty. I changed jobs. Still the same, I was dreading in boredom and rat race. Needing to "impress" people or feeling that I need to be there for work 9-5 each day without feeling motivated was awful. I quit after a few months willing to take some time off. What am I doing? I have no purpose and my day to day means nothing. I like seeing my friends and family, traveling, etc but the regular life has a lot more things and time to fill. I can't make this for 30 more years.

Fortunately, tech and ML pays high salaries so I felt safe. I could take some time off, I can always look for other job, etc. But the feeling of "I need to do something different" was there. I started looking for other opportunities, looking for purpose in jobs, looking for motivation. I was in despair as I understood that I needed to do my own way and no one was doing it for me.

Ironically, I took another job quit quickly because it paid much more and I could just try it, maybe I was just in the wrong place. Seemed to work better, it was less of a charade, I had more time for crafting and working, I got some motivation back as I felt at least more comfortable daily. Not that I had purpose really but working was at least with less pressure. The time zone differences helped because I didn't feel like 9 hours a day I was supposed to be there, with someone 1 click away of sending a slack message to me. I stayed two years. Earning good money compared to my spending, saving and "happy". Eventually everything started to decay, the job was less free, the business request were more urgent and I quickly lost motivation and started to fall in despair again. What am I doing? Who cares about this parquet file with funny numbers?

I quit, this time for good. Despair again, realization, I need to find something for myself, or try at least. I fear for its complexity, I fear if I can do it, if I can make a living in another way. I'm taking some time to rest but I know eventually I will need to figure things out. On my own, making my own way, and I see how no one will do it for me.

As a side note, finding purpose, finding your way, etc was also a thing when I thought about romantic relationships. Finding someone that you really care about and that cares about you is not easy at all and it doesn't happen magically as I thought as a kid. I am not going to expand about it here but it was another topic that made me realize about our own fate and efforts.

Side side note, retirement is another one. I don't see my national retirement plan as something you can live off. I see my family, that luckily does well but despite that nothing is granted and we can support the elder in my family because our own safety net. A younger me didn't know that but as my twenties went by I started to see a lot more of adult/real life. Eye opening period of my life.

How to be lucky

"Okay, I'm focused on getting X, but let's not forget to read the headlines."

High Agency and be the plumber

Reminder to myself, focus on bringing solutions and not the shiny new tool. I think I'm good at high agency or at being responsible and wanting to finish what I commit to.
I need to work on focusing on which solution am I bringing. Finding the right niche and actually bringing value. In my case with sportsjobs.online, not just throwing things into it, but making it clear what I am providing. For other topics related to ML and LLM I need to decide what I'll focus on this time coming in.

Impostor Syndrome I'm the classical example of not trusting myself and thinking stuff like this.

but at the end of the day, you must just think I have shit taste and that you've somehow tricked me into thinking you're good when you're an impostor? Right?

I need to stop with that. Quitting now I had plenty of nice words for every colleague, including the desire from managers that I stay or I come back if I get bored of not having a job. Of course, as I write that, I think in the back of my mind of exceptions and that technical colleagues were less effusive as my managers and etc. but all of that is probably not true. I'm going against all evidence and just guessing and making things in my head.

How to Be Good at Many Things Consistency. Every one says this. I need to do it the right way now that I'll have the time. No excuses.

Do things, practice, keep going for it. Everything will be easier and quicker.

And I need to be grateful for what I have and how lucky I am to be able to get a sabbatical time to try to change the path I'm going.

Weighted regression

Weighted regression consists on assigning different weights to each observation and hence more or less importance at the time of fitting the regression.

On way to look at it is to think as solving the regression problem minimizing Weighted Mean Squared Error(WSME) instead of Mean Squared Error(MSE)

\(\(WMSE(\beta, w) = \frac{1}{N} \sum_{i=1}^n w_i(y_i - \overrightarrow {x_i} \beta)^2\)\) Intuitively, we are looking fot the coefficients that minimize MSE but putting different weights to each observation. OLS is a particular case where all the \(w_i = 1\)

Why doing this? A few reasons (Shalizi 2015. Chapter 7.1)

  • Focusing Accuracy: We want to predict specially well some particular points or region of points, maybe because that's the focus for production or maybe because being wrong at those observations has a huge cost, etc. Using Weighted regression will do an extra effort to match that data.

  • Discount imprecision: OLS returns the maximum likelihood estimates when the residuals are independent, normal with mean 0 and with constant variance. When we face non constant variance OLS no longer returns the MLE. The logic behind using weighted regression is that makes no sense to pay equal attention to all the observations since some of them have higher variance and are less indicative of the conditional mean. We should put more emphasis on the regions of lower variance, predict it well and "expect to fit poorly where the noise is big".
    The weights that will return MLE are \(\frac{1}{\sigma_i^2}\)

  • Sampling bias: If we think or know that the observations in our data are not completely random and some subset of the population might be under-represented (in a survey for example or because of data availability) it might make sense to weight observation inversely to the probability of being included in the sample. Under-represented observations will get more weights and over-represented less weight.
    Another similar situation is related to covariate shift. If the distribution of variable x changes over time we can use a weight designed as the ratio of the probability density functions.

    "If the old pdf was p(x) and the new one is q(x), the weight we'd want to is \(w_i=q(x_i)/p(x_i)\)

  • Other: Related to GLM, when the conditional mean is a non linear function of a linear predictor. (Not further explained in the book at this point)

Is there a scenario where OLS is better than Weighted regression? Assuming we can compute the weights.

Example.

First we will see the impact of using weighted regression, using a simulated scenario where we actually know the variance of the error of each observation. This is not realistic but useful to see it in action.

library(tidyverse)

We generate 1000 datapoints with a linear relation between y and x. Intercept = 0, slope = 5. We let the variance of the error depend on the value of x. Higher values of x are associated with higher values of the variance of the error.

set.seed(23)
n=1000
x = runif(n,0,10)
error = rnorm(n,0, x/1.5)
df = data.frame(x)
df = df %>% mutate(y = 5*x + error)
Visually..

ggplot(data=df, aes(x=x, y=y)) +
  geom_point(alpha=0.3) + 
  # geom_smooth(color="blue") +
  # geom_smooth(method = "lm", mapping = aes(weight = (1/sqrt(x)^2)),
  #             color = "red", show.legend = FALSE) +
  NULL

Image

Linear regression
ols = lm(formula = "y~x", data=df)
summary(ols)


## 
## Call:
## lm(formula = "y~x", data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.868  -1.720  -0.137   1.918  14.722 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.19192    0.24278   0.791    0.429    
## x            4.95585    0.04148 119.489   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.855 on 998 degrees of freedom
## Multiple R-squared:  0.9347, Adjusted R-squared:  0.9346 
## F-statistic: 1.428e+04 on 1 and 998 DF,  p-value: < 2.2e-16
We get an intercept of 0.19, non-significant when the actual value is 0 and a slope of 4.96 when the actual value is 5.

Weighted linear regression
wols = lm(formula = "y~x", data=df, weights = (1/sqrt(x)^2) )
summary(wols)

## 
## Call:
## lm(formula = "y~x", data = df, weights = (1/sqrt(x)^2))
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8880 -0.8601 -0.0016  0.8936  4.6535 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.001483   0.030072   0.049    0.961    
## x           4.993473   0.021874 228.286   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.498 on 998 degrees of freedom
## Multiple R-squared:  0.9812, Adjusted R-squared:  0.9812 
## F-statistic: 5.211e+04 on 1 and 998 DF,  p-value: < 2.2e-16

We get an intercept of 0, non-significant too but much closer to 0 and with lower standard error and a slope of 4.99 also much closer to the actual value of 5 and with lower standard error.

Conclusion: if we know the right weights we can get better estimates from a linear regression in case of heteroscedasticity.

Inference is not valid in the dataset used for model selection.

Let's say we have a dataset and we want to fit a model to it and do some inference such as obtaining the coefficients and look for their confidence intervals.

For such a task we would first need to find a model that we think approximates to the real data generating process behind the phenomenon.
This will be the model selection step.
Then we would look at the output of our model and get the standard error of the coefficients or calculate the confidence interval or any other similar task. This will be the inference step.

The issue here is that, if we don't know the true model and we do model selection, our own model will be a random object. Why? Because the particular dataset we are using is also a set of random variables. Other datasets might return another model formula as the best between our options since that particular dataset would have other observations and particularities.

Main problem:

since we are selecting a model based on a particular dataset, the standard errors and p-values will be smaller than then actual ones.

"That means there is some extra randomness in your estimated parameters (and everything else), which isn't accounted for by formulas which assume a fixed model.
This is not just a problem with formal model-selection devices like cross-validation. If you do an initial, exploratory data analysis before deciding which model to use - and that's generally a good idea - you are, yourself, acting as a noisy, complicated model-selection device" (Sharizi 2017)

The most straightforward way to deal with this (if you are using independent observations) is to split the data, do model selection in one part and then fit the best model in the other part. Your second fit will be the one useful for inference.
You could fit the model to the full data but that would include the part used for model selection and you would still get false, overconfident standard errors.

Let's see an example.
We will generate data following a somewhat "complicated" model with interactions. We will split the data in two equal size parts. One for model selection and one for inference.
We will then fit a couple formulas to model selection part and pick the one with the minimum RMSE. We will compare the standard errors obtained in the model selection part and the ones obtained fitting that model to the inference part.

Thanks to BrunoRodrigues for this post that I used as guideline to fit models with Cross Validation in R.

We start by generating the data, including interactions.

set.seed(1)
N = 5000
b0 = 4
b1 = 1
b2 = 2
b3 = 3
b4 = 4
b5 = 5

x1 = runif(N, 0, 10)
x2 = rnorm(N, 20, 3)
x3 = runif(N, 20, 40)
error = rnorm(N, 0, 200)

y = b0 + b1*x1 + b2*x2 + b3*x3 + b4*x1*x2 + b5*x2*x3 + error


df = tibble(y, x1, x2, x3)

We do the first split, df_selection will be the one used to try different models and pick one.
df_inference will be used to do the actual inference given the model selected.

prop = 0.5

selection_inference_split = initial_split(df, prop=prop)

df_selection = training(selection_inference_split)
df_inference = testing(selection_inference_split)

To select a model using df_selection we will use Cross validation to try to get the model that best generalizes.
We will generate 30 split of 70% of the data and use the other 30% to calculate RMSE metric.

validation_data <- mc_cv(df_selection, prop = 0.7, times = 30)

We create two functions, my_lm() will run a linear regression for the training part of each split of CV and get the prediction for the testing part of each split. We will run this for a couple of formulas.
return_model will fit the model to the whole training data to extract the parameters and standard errors we get if we use the same dataset that was used to do model selection.

my_lm <- function(formula, split, id){

    analysis_set = analysis(split)  

    model <- lm(formula = formula, data=analysis_set)

    assessment_set <- assessment(split)


    pred = tibble::tibble("id" = id,
        "formula" = formula,
        "truth" = assessment_set$y,
        "prediction" = unlist(predict(model, newdata = assessment_set)))

}


return_model = function(formula){


    model <- lm(formula = formula, data=df_selection)
    output = list(model$coefficients, summary(model))

}

We will try 5 formulas. The first one is the actual data generating process and should the best in terms of RMSE. We will exclude that one for model selection since the aim of this is to simulate a scenario where we don't know the actual formula behind the data. We calculate it just for reference but we will pick one of the other 4 models for inference.

formulas = list("y ~ x1 + x2 +x3 + x1*x2 + x2*x3", 
                "y ~ .", 
                "y ~ x1 + x2", 
                "y ~ x1 + x2 + x3 + x1*x2",
                "y ~ x1 + x2 + x3 + x2*x3")
results = data.frame()

models = list()
for (formula in formulas){

results_selection <- map2_df(.x = validation_data$splits,
                           .y = validation_data$id,
                           ~my_lm(formula = formula, split = .x, id = .y))

model = return_model(formula)


results = rbind.data.frame(results, results_selection)
models = c(models, model )

}

We retrieve the mean RMSE across the splits, calculated in the test part of each split.
We can see that the real model is the best in terms of RMSE. Between the others, we can see that the one including the x2:x3 interaction is the best. So, we will keep that one as our "model selected"

results %>%
    group_by(id, formula) %>%
    rmse(truth, prediction) %>%
    group_by(formula) %>%
    summarise(mean_rmse = mean(.estimate)) %>%
    as.data.frame()

##                           formula mean_rmse
## 1                           y ~ .  219.4756
## 2                     y ~ x1 + x2  625.0173
## 3        y ~ x1 + x2 + x3 + x1*x2  217.3185
## 4        y ~ x1 + x2 + x3 + x2*x3  198.9802
## 5 y ~ x1 + x2 +x3 + x1*x2 + x2*x3  196.4747

We can check the parameters and the standard errors when fitted to the whole selection dataset.

## 
## Call:
## lm(formula = formula, data = df_selection)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -854.07 -132.91   -0.97  137.65  714.53 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -245.1084   138.9684  -1.764   0.0779 .  
## x1            82.7919     1.3540  61.148   <2e-16 ***
## x2            15.7118     6.8628   2.289   0.0221 *  
## x3            -1.5177     4.5626  -0.333   0.7394    
## x2:x3          5.1676     0.2256  22.906   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 199.1 on 2495 degrees of freedom
## Multiple R-squared:  0.947,  Adjusted R-squared:  0.9469 
## F-statistic: 1.115e+04 on 4 and 2495 DF,  p-value: < 2.2e-16

And let's see what happens if we fit the same model to the inference set.

model_test = lm(formula=formulas[[5]], data=df_inference)

summary(model_test)


## 
## Call:
## lm(formula = formulas[[5]], data = df_inference)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -656.47 -138.67   -5.64  130.21  773.99 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -438.7059   140.2618  -3.128 0.001782 ** 
## x1            81.4724     1.3622  59.812  < 2e-16 ***
## x2            23.4856     6.9475   3.380 0.000735 ***
## x3             3.6309     4.5942   0.790 0.429417    
## x2:x3          4.9750     0.2275  21.869  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 199.4 on 2495 degrees of freedom
## Multiple R-squared:  0.9473, Adjusted R-squared:  0.9472 
## F-statistic: 1.121e+04 on 4 and 2495 DF,  p-value: < 2.2e-16
First we can see that the parameters have changed a bit.
In second place we can see that the standard errors are generally bigger in comparison to the parameter for the inference set and will generate a wider confidence interval.

ggplot(data = ratio_df) +
  geom_point(aes(x=parameter, y=ratio, col=set), size=3) +
  theme(legend.title = element_blank()) +
  theme_light() +
  xlab("") +
  ylab("Ratio") +
  ggtitle("Absolute ratio between SD and Estimate")

Image

My idea is to add a plot with the confidence intervals so the effect can be seen directly but I don't have the time today. Anyways, it is clear that the standad error to parameter ratio is bigger in the inference set, showing that the inference in the same dataset as model selection is invalid as it is overconfident in the results.

Remarks on R2

R2 depends on the variance on the variance of the predictors

Quoting from Shalizi[^1] Assuming a true linear model
$$ Y = aX + \epsilon$$
and assuming we know \(a\) exactly.
The variance of Y will be \(a^2\mathbb{V}[X] + \mathbb{V}[\epsilon]\).
So \(R^2 = \frac{a^2\mathbb{V}[X]}{a^2\mathbb{V}[X] + \mathbb{V}[\epsilon]}\)
This goes to 0 as \(\mathbb{V}[X] \rightarrow 0\) and it goes to 1 as \(\mathbb{V}[X] \rightarrow \infty\). "It thus has little to do with the quality of the fit, and a lot to do with how spread out the predictor variable is. Notice also how easy it is to get a high \(R^2\) even when the true model is not linear!"

Below a quick comparison between two linear relationships, one with much higher variance than the other in the predictor.
Added a different constant for better display in plot.

library(tidyverse)

x1 = rnorm(1000, mean=0, sd=1)
x2 = rnorm(1000, mean=0, sd=10)
error = rnorm(1000, mean=0, sd=0.5)

y1 = x1 + error
y2 = 10 + x2 +  error

df = data.frame(x1,x2,y1,y2)

model1 = lm("y1 ~ x1")

## Error in eval(predvars, data, env): object 'y1' not found

model2 =  lm("y2 ~ x2")

## Error in eval(predvars, data, env): object 'y2' not found

Linear Smoothers

Linear regression as smoothing

Let's assume the DGP (data generating process) is: $$ Y = \mu(x) + \epsilon$$ where \(\mu(x)\) is the mean Y value for that particular x and \(\epsilon\) is an error with mean 0.

When running OLS we are trying to approximate \(\mu(x)\) with a linear function of the form \(\alpha + \beta x\) and trying to retrieve the best \(\alpha\) and \(\beta\) minimizing the mean-squared error.

The conclusions don't change but the math gets easier if we assume both X and Y are centered (mean=0).
With that in mind we can write down the MSE and optimize to get the best parameters.

\[MSE(\alpha, \beta) = \mathbb{E}[(Y - \alpha - \beta X)^2] \\ = \mathbb{E}[\mathbb{E}[(Y - \alpha - \beta X)^2 | X]] \\ = \mathbb{E}[\mathbb{V}[Y|X]] + \mathbb{E}[Y- \alpha - \beta X | X])^2] \\ = \mathbb{E}[\mathbb{V}[Y|X]] + \mathbb{E}[(\mathbb{E}[Y- \alpha - \beta X | X])^2]\]

Deriving with respect to \(\alpha\) and \(\beta\) for optimization..
The first term can be dropped since doesn't include any parameter.

$$\frac{\partial MSE}{\partial \alpha} = \mathbb{E}[2(Y - \alpha - \beta X)(-1)] \ \mathbb{E}[Y - a - b X] = 0 \ a = \mathbb{E}[Y] - b \mathbb{E}[X] = 0 $$ when Y and X are centered..

and $$\frac{\partial MSE}{\partial \beta} = \mathbb{E}[2(Y - \alpha - \beta X)(-X)] \ \mathbb{E}[XY] - b\mathbb{E}[X^2] = 0 \ b = \frac{Cov[X,Y]}{\mathbb{V}[X]} $$

The optimal beta is a function of the covariance between Y and X, and the variance of X.

Putting together \(a\) and \(b\) we get \(\mu(x) = x \frac{Cov[X,Y]}{\mathbb{V}[X]}\)

Replacing with the values from the sampled data we get an estimation of \(a\) and \(b\).

Remember they are 0 centered so variance and covariance get simplified.

\[ \hat a = 0 \\ \hat b = \frac{\sum_i y_i x_i}{\sum_i x_i^2}\]

With all this we can see how OLS is a smoothing of the data.
Writing in terms of the data points:
$$\hat \mu(x) = \hat b x \ = x \frac{\sum_i y_i x_i}{\sum_i x_i^2} \ = \sum_i y_i \frac{x_i}{\sum_j x_j^2} x \ = \sum_i y_i \frac{x_i}{n \hat \sigma_x^2} x $$ where \(\hat \sigma_x^2\) is the sample variance of X.
In words, our prediction is a weighted average of the observed values \(y_i\) of the dependent variable, where the weights are proportional to how far \(x_i\) is from the center (relative to the variance), and proportional to the magnitude of \(x\). If \(x_i\) is on the same side of the center as \(x\), it gets a positive weight, and if it's on the opposite side it gets a negative weight. (Shalizi 2017)

If \(\mu(x)\) is really a straight line, this is fine, but when it's not, that the weights are proportional to how far they are to the center and not the point to predict can lead to awful predictions.

Alternative smoothers

For that, other methods smooth the data in another ways to help mitigate that.

As quick examples, we have KNN regression where the smoothing is done using only close observations to the one to predict (and getting quite noisy since depend a lot on the sample points around a small area).

Kernel smoothers are a variant where depending on the kernel selected we get different smoothing. The main idea is that we use a windowd of data with the idea of putting more weight to points close to the one to predict. Could be Gaussian weight around X for example, or uniform around a window. Note this is different than KNN regression since we do not take the average of those points, we get a regression for that area.
A nice thing about this smoothers (and KNN regression) is that if we want to predict points far from the training data we won't get a linear extrapolation as with OLS but it will be pushed towards the closest data points we had in training.

Bias Variance Tradeoff

Mean squared error (MSE) is a measure of how far our prediction is from the true values of the dependent variable. It's the expectation of the squared error.

The squared error being:

\[(Y - \hat \mu(x))^2\]

where Y is the true value and \(\hat \mu(x)\) is the prediction for a given x.

We can decompose it into:

\[(Y - \hat \mu(x))^2 \\ = (Y - \mu(x) + \mu(x) - \hat \mu(x)^2) \\ = (Y - \mu(x))^2 + 2(Y - \mu(x))(\mu(x) - \hat \mu(x)) + (\mu(x) - \hat \mu(x))^2\]

So, that's the squared error. The MSE is the expectation of that.

The expectation is a linear operator so we can apply it independently to different terms of a summation.
The expectation of the first term is the variance of the error intrinsic to the DGP.
The second term goes to 0 because involves \(E(Y-\mu(x))\) that is the expectation of the error and that's equal to 0.
The third term reamins as it is since doesn't involve random variables.

\[MSE(\hat \mu(x)) = \sigma^2_x + (\mu(x) - \hat \mu(x))^2\]

This is our first bias-variance decomposition. The first term is the intrinsic difficulty of the problem to model, is the variance of the error and can not be reduced, it is what it is.
The second term is how off our predictions are regarding the true expected value for that particular X.

This would be fine if we wouldn't need to consider \(\hat \mu(x)\) a random variable itself, since it is dependent on the specific dataset we are using. Given another dataset our estimation would be different despite using the same model methodology.
What we actually want is the MSE of the method used \(\hat M\) and not only the result of a particular realization.

\[MSE(\hat M_n(x)) = E[(Y - \hat M_n(X))^2 | X=x] \\ = ... \\ = \sigma^2_x + (\mu(x) - E[\hat M_n(x)])^2 - V[\hat M_n(x)] \]

This is our 2nd bias-variance decomposition.
The first term is still the irreducible error.
The second term is the bias of using \(\hat M_n\) to approximate \(\mu(x)\). Is the approximation bias/error.
The third term is the variance of the estimate of the regression function. If our estimates have high variance we can have large errors despite using an unbiased approximation.

Flexible methods will be able to approximate \(\mu(x)\) closely, however usually using more flexible methods involve increasing the variance of the estimate. That's the bias-variance tradeoff. We need to evaluate how to balance that, sometimes including some bias reduce much more the error by decreasing the variance.
Usually larger N decreases the MSE since it decreases bias and variance error.

Reference

Based on 1.4.1 from Advanced data analysis from a elementary point of view.

Spark and Pyspark

What's Spark?

prueba The definition says:

Spark is a fast and general processing engine compatible with Hadoop data. It can run in Hadoop clusters >through YARN or Spark's standalone mode, and it can process data in HDFS, HBase, Cassandra, Hive, and any >Hadoop InputFormat. It is designed to perform both batch processing (similar to MapReduce) and new >workloads like streaming, interactive queries, and machine learning.

Basically is a framework to work with big amounts of data stored in distributed systems instead of just one machine. This allows parallelization and hence much faster calculations.
It's biggest difference with plain Hadoop is that Spark uses RAM to process data while Hadoop doesn't.

Not being a data engineer myself I can tell you that you can use Spark to work with data stored in HDFS, S3 buckets or a data lake for example. All distributed systems.

Since those usually store huge big amount of data you can see how all this relate. The use case I have been exposed to, as a data scientist, is to query this distributed data and process it before using it for some purpose (modeling, reporting, etc).

How to use it?

I haven't deployed a distributed storage system myself but I think it's safe to assume that amount of data is gathered in big organizations and probably some data engineer has already done all the setup. You just want to access the data from an environment connected to the spark cluster.

There are several languages that can interact with Spark. Scala is the original one but you could use Java or Python. As data scientist we are probably more familiar with Python so I will show you Pyspark

Pyspark

Pyspark is an API to work with Spark using Python. In order to run you need also Java installed and Apache Spark. In our fictional organization a data engineer might have set up a server with Jupyter notebooks linked to the data lake and with all the dependencies.

There are probably ways to connect to the remote spark server from your local machine but I haven't done that.

So, Pyspark allows you to query the datalake/bigdata storage from a jupyter notebook and then convert that to a Pandas Dataframe and work as you are used to.

Spark/Pyspark has a particular syntax that is quite clear but has some particularities based on the parallelization notion. For example, many functions don't actually retrieve all the data, that only happens when you decide to. For example show() or collect() do retrieve the data (and can take a while if you are working with a lot of data) while filter() or withColumn() don't.

Another thing to notice is that you will need to create/initiate a sparkContext before actually being able to query data.

To understand this and have a good amount of examples regarding the functions and syntax I highly recommend THIS SITE.

How to practice?

You can practice Pyspark queries and scripts by installing Pyspark in your local machine despite not having a cluster running distributed data. With Pyspark installed you can create some data and use it as it was real.
You will be able to use all the functions and check them by yourself.

How to install it? You can check THIS GUIDE FOR WINDOWS.

I have struggled a bit to make it work so these are some things I learned during the way.

  • I have downloaded Java 8 since that's what the guide says and use that at my current organization.
  • To avoid creating an account in Oracle to download Java you can check THIS SOLUTION.
  • When creating Environment variables avoid blank spaces
  • If Pyspark doesn't run because can't find Java. Check the %JAVA_HOME% path.
  • If the error is related to missing Python3 , check the %PYTHONPATH% and create in the anaconda path a copy of python.exe but rename it python3.exe