Criando um dataset:

set.seed(150)
data <- data.frame(x = rnorm(50, mean = 50, sd = 10),
                   randomico = sample(c(-10:10), 50, replace = TRUE))
data$y <- data$x + data$randomico
head(data, 10)

Calculando a correlação entre X e Y.

correlacao <- cor(data$x, data$y, method = 'pearson')
correlacao
## [1] 0.8869656

Interpretando o coeficiente:
r=0; Não há correlação entre as variáveis.
r=+1; Correlação positiva perfeita entre as variáveis.
r=-1; Correlação negativa perfeita entre as variáveis.
r= 0 to 0.30; correlação pouco significativa entre as variáveis.
r=0.30 to 0.50; correlação moderada entre as variváeis.
r=0.50 to 1 correlação forte entre variáveis.

set.seed(150)
valoresX <- rnorm(50, mean = 50, sd = 10)
randomico <- sample(c(10:30), 50, replace = TRUE)
data <- data.frame(x = rep(valoresX, 2),
                   randomico = rep(randomico, 2),
                   Categoria = rep(c("Um","Dois"), each = 50))
data$y[data$Categoria=="Um"] <- 20 + data$x[data$Categoria=="Um"]/data$randomico[data$Categoria=="Um"]
data$y[data$Categoria=="Dois"] <- 20 + data$x[data$Categoria=="Dois"]/(5*data$randomico[data$Categoria=="Dois"])
data
correlacaoUm <- cor(data$x[data$Categoria=="Um"], data$y[data$Categoria=="Um"], method = 'pearson')
correlacaoDois <- cor(data$x[data$Categoria=="Dois"], data$y[data$Categoria=="Dois"], method = 'pearson')

Normalmente se comete o equívoco de associar a correlação de Pearson com a inclinação (slope) de uma linha entre as duas variáveis, como no caso do coeficiente angular de uma regressão linear.
No entanto, o coeficiente de Pearson mede apenas a força da relação entre as variãveis em consideração. Isso pode ser demonstrado abaixo ao verificar os coeficientes obtidos e as retas de uma regressão linear aplicada ao mesmo conjunto de dados.

correlacaoUm
## [1] 0.5476772
correlacaoDois
## [1] 0.5476772
library(ggplot2)
gg <- ggplot(data, aes(x, y, colour = Categoria))
gg <- gg + geom_point()
gg <- gg + geom_smooth(alpha=0.3, method="lm")
gg
## `geom_smooth()` using formula 'y ~ x'

Categoria Um:

fit = lm(x~y, data = data[data$Categoria == "Um", ])
summary(fit)
## 
## Call:
## lm(formula = x ~ y, data = data[data$Categoria == "Um", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.180  -5.479  -1.249   5.951  20.364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -67.710     26.164  -2.588   0.0127 *  
## y              5.216      1.150   4.535 3.85e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.224 on 48 degrees of freedom
## Multiple R-squared:    0.3,  Adjusted R-squared:  0.2854 
## F-statistic: 20.57 on 1 and 48 DF,  p-value: 3.848e-05

Categoria Dois:

fit = lm(x~y, data = data[data$Categoria == "Dois", ])
summary(fit)
## 
## Call:
## lm(formula = x ~ y, data = data[data$Categoria == "Dois", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.180  -5.479  -1.249   5.951  20.364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -484.996    118.152  -4.105 0.000156 ***
## y             26.080      5.751   4.535 3.85e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.224 on 48 degrees of freedom
## Multiple R-squared:    0.3,  Adjusted R-squared:  0.2854 
## F-statistic: 20.57 on 1 and 48 DF,  p-value: 3.848e-05