
Above is the scatter plot of my data. The spearman's rho was 0.66 which was significant (p<0.005). I draw the fitting line using polynomials. My question: is it OK to choose arbitrary function for fitting line in this case? If so, what would be the best function for the curve fitting?

Above was drawn according to the below recommendation by Glen_b (Thanks a lot). Below is the R code for it.
install.packages('uniReg')
library(uniReg)
x <- read.delim('clipboard', header=FALSE)
y <- read.delim('clipboard', header=FALSE)
fit <- unireg(x$V1,y$V1,g=5,sigma=2)
z=seq(min(x$V1),max(x$V1),length=201)
uf=with(x,fit)
plot(x$V1, y$V1,type = 'p', col = 'blue', pch=16, xlab="LLG", ylab="TSH")
lines(z,uf$unimod.func(z),type = 'l', col = 'red' )
uniregfunction fits both (among others) -- the argumentconstr="isotonic"should do that, but I'm not sure it's better than your earlier fit; you may also need to changepenaltyand some of the other arguments to get a suitable fit (the default looks too "stiff"). If you only need this one fit, there's no problem; the fit you have is isotonic. I should probably make my comments into some form of answer. – Glen_b Jun 30 '15 at 10:19