How to analyze relationship between data pairs Please change the title if it's not accurate, I might be lacking terminology.
I have two variables ($x,y$) which are directly correlated by pair. For example I have sampled 50 trees in my study area and for each tree I have tree height and annual precipitation at that location. This is what I mean by they are organized by pair. So I plot $x$ vs $y$ and I would like to know if there is any correlation between the two; does annual precipitation predict tree height, or the other way around?
I have used R to test if there is any linear correlation with the function lm(y~x) and have $p<.05$ and $R^2=0.6$, but some of my data seem as if they might be better fitted with a non-linear curve. What are some other types of statistical analysis I should use for such a problem?
 A: Kernel regression and smoothing splines will both fit non-linear regression functions very nicely. There's a longer list of methods, with examples, in this R-bloggers post by Wesley. 
A: To visualize nonlinear data, try: 
library(ggplot2)
ggplot(rndf, aes(x=vnum1, y=vnum2, group=vbin1, color=vbin1))+stat_smooth()

Where vnum1 and vnum2 are your x and y and vbin1 is the species or other grouping column. 

Different options can be set in stat_smooth() function. 
The shaded areas indicate 95% CI. Overlapping areas would indicate non-significant difference. 
A: Rank correlation coefficients such as Kendall's or Spearman's can be useful to measure the strength of non-linear monotonic relationships, i.e. those that are increasing or decreasing, but not e.g. U-shaped. Calculating confidence intervals for their population values is often considered more useful than simply testing whether they're equal to zero.
If you want to predict tree height from precipitation (or vice versa) there are many approaches (see @Lizzie's links). The most basic extension of simple linear regression is simply to represent the predictor $x$ as a polynomial & estimate a coefficient $\beta_k$ for each term $k$:
$$\operatorname{E}Y =\beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + ...$$
where $\operatorname{E}Y$ is the predicted response. With just one predictor you might use F-tests to see how many terms are necessary (too many compared to the number of observations & you'll over-fit the model to the sample data).
