Find "best" range of one value based on correlated value Scenario: I want to give advice on how long a post should be for maximum number of interactions (lets call "Likes"). I have Description Length and Likes for each post. 
(I've excluded posts with no Likes, is this valid or a statistical faux pas?)
Scatter plot of Likes v Description Length:

I could just look at the plot and say 800–3500 is the best length range, with the optimum being around 2000 characters. But how do I find that out statistically?
So I guess I should:


*

*Exclude outliers / smooth data somehow (use bins?)

*Find the "best" number of Likes? Or best Lengths for the most Likes?

*Find the range of Lengths that falls within 1 standard deviation of that (how)?

*Or maybe I want to fit a curve and find the highest point?


I'm sorry this is probably a pretty basic problem, I'm just not really sure what to search for ("find corelation between data sets within one standard deviation of best value" didn't help). So names of techniques and links to articles might be enough. Thanks!
Data set: https://cl.ly/lPro
 A: First: You should definitely include the posts with zero likes in the data. They are informative, and without them the results could be biased. That being said, I will show one possible approach using the data you have posted. 
The count of likes, as a count variable, could be taken as poisson distributed. However, with such large counts and many possibilities for unobserved heterogeneity, I will use a generalized linear model (glm) with a quasipoisson family. Then I will use a spline model, plotting the predictions from that model (with red), and predictions from some bootstrap resamples to show the uncertainty in the predictions, informally. The result is:
 
The predictions from the bootstrap in orange. It is quite clear from this plot that not much can be said about the optimal post length. 
Below some code in R:
library(splines)
set.seed(7*11*13) # my public seed 
mod  <-  glm(likes ~ ns(length, df=5), data=LLdat, family=quasipoisson())  

pred  <-  predict(mod, type="response")   

LLdat  <-  within(LLdat, pred <- pred)

o <- with(LLdat, order(length))
LLdat.sort  <-  with(LLdat, data.frame(length=length[o], pred=pred[o], likes=likes[o]))   

### Bootstrapping
### We do it by hand:

R <- 15
n <- 743
o <- with(LLdat, order(length))
newdata <- LLdat[o, ]
m <- NROW(newdata)
predictions <-  matrix(0.0, m, R)

with(LLdat.sort, plot(length, likes, log="y", col=scales::alpha("blue",0.4), main="'likes' versus post length", cex.main=sqrt(2)))  
with(LLdat.sort, lines(length, pred, col="red"))   
for (r in 1:R) { ind <- sample(1:n,n,replace=TRUE) ; mod <- glm(likes ~ ns(length, df=5), data=LLdat[ind,], family=quasipoisson()) ; predictions[,r] <-  predict(mod, type="response",newdata=newdata) }
for (r in 1:R) {lines(newdata[,1], predictions[,r], col="orange") } 

