library(corrplot)
library(reshape2)
library(ggplot2)
library(caret)
## Loading required package: lattice

Os dados conciste de vol: volume do câncer weight: peso do paciente age: idade do paciente bph: hiperplasia prostática benigna svi: invasão das vesículas seminais cp: penetração capsular gleason: escore Gleason pgg45: percentagem escore Gleason 4 ou 5 psa: antígeno específico da próstata (esta é a variável resposta).

prostate.data <- read.delim("~/Dropbox/UFCG/AD2/M3/prostate.data.txt")

No gráfico abaixo é mostrado as distribuições das variáveis em que no eixo x representa o valor correspondente da variável e no eixo y a quantidade daquele determinado valor que apareceu nos dados

df <- melt(prostate.data)
## Using train as id variables
ggplot(df,aes(x = value)) + 
  facet_wrap(~variable, scales = "free_x") + 
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Como é possível observar algumas variáveis apresentam uma distribuição normal dos dados, ou seja existe um equilíbrio entre os menores valores e os maiores valores. Algumas variáveis tem um intervalo de valores muito pequeno como é o caso das variáveis gleason e svi.

Primeiramente eu pensei em realizar uma normalização dos dados, mas acredito que tal processo não melhoraria os resultados de forma significativa.

No gráfico abaixo podemos ver a correlação entre diversas variáveis, em outras palavras o qual similar são tais varíaveis. Por exemplo, a variável gleason tem um comportamento muito parecido com a variável pgg45.

M <- cor(prostate.data)
corrplot(M,type="lower")

Se elas têm comportamento parecido então só precisamos de apenas uma das duas.

O mesmo ocorre com as variáveis lcp e lcavol, e lcp e svi.

A pergunta que fica é: Qual das duas devo eliminar? Eu optei pela variável que tem menor correlação com a variável resposta, ou seja a variável que queremos prever.

Como podemos ver nos gráficos seguintes, há uma comparação com todos os dados e após a remoção das varíaveis lcp e gleason Abaixo temos o sumário de uma regressão linear que de forma geral é um método estatístico para prever valores de uma varíavel a partir do processamento de outras variáveis.

train <- prostate.data[prostate.data$train==TRUE,]
test <- prostate.data[prostate.data$train==FALSE,]
modeloA <- lm(lpsa ~ lcavol
              + lweight
              + age
              + lbph
              + svi
              + pgg45
              + lcp
              + gleason,data = train)
summary(modeloA)
## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + pgg45 + 
##     lcp + gleason, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.64870 -0.34147 -0.05424  0.44941  1.48675 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.429170   1.553588   0.276  0.78334    
## lcavol       0.576543   0.107438   5.366 1.47e-06 ***
## lweight      0.614020   0.223216   2.751  0.00792 ** 
## age         -0.019001   0.013612  -1.396  0.16806    
## lbph         0.144848   0.070457   2.056  0.04431 *  
## svi          0.737209   0.298555   2.469  0.01651 *  
## pgg45        0.009465   0.005447   1.738  0.08755 .  
## lcp         -0.206324   0.110516  -1.867  0.06697 .  
## gleason     -0.029503   0.201136  -0.147  0.88389    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7123 on 58 degrees of freedom
## Multiple R-squared:  0.6944, Adjusted R-squared:  0.6522 
## F-statistic: 16.47 on 8 and 58 DF,  p-value: 2.042e-12
predictions <- predict(modeloA,test)
axisRange = extendrange(c(test$lpsa,predictions))
plot(test$lpsa,predictions,main="Todas as Variaveis")
abline(0,1,col="blue",lty=2,lwd=2)

Abaixo podemos observar o erro residual, ou seja o qual as predições se afastam do valor real.

lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)
##     RMSE Rsquared 
##    0.722    0.505
modeloB <- lm(lpsa ~ lcavol
              + lweight
              + age
              + lbph
              + svi
              + pgg45,data = train)
summary(modeloB)
## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + pgg45, 
##     data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.79922 -0.30470 -0.05402  0.41301  1.45256 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.428057   1.042469   0.411  0.68281    
## lcavol       0.488440   0.096658   5.053 4.35e-06 ***
## lweight      0.613739   0.223090   2.751  0.00784 ** 
## age         -0.017273   0.013323  -1.296  0.19979    
## lbph         0.154229   0.071066   2.170  0.03396 *  
## svi          0.553158   0.282538   1.958  0.05491 .  
## pgg45        0.005358   0.003702   1.447  0.15306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7211 on 60 degrees of freedom
## Multiple R-squared:  0.676,  Adjusted R-squared:  0.6436 
## F-statistic: 20.86 on 6 and 60 DF,  p-value: 4.845e-13
predictions <- predict(modeloB,test)
axisRange = extendrange(c(test$lpsa,predictions))
plot(test$lpsa,predictions,main="Apos exclusao de gleason e lcp")
abline(0,1,col="blue",lty=2,lwd=2)

Podemos observar que o erro residual diminuição em relação ao nosso modelo anterior.

lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)
##     RMSE Rsquared 
##    0.678    0.561

Utilizando uma abordagem de retirada de variáveis chegamos ao modelo final a qual temos um erro residual inferior aos dois modelos anteriores.

modeloC <- lm(lpsa ~ lcavol
              + lweight
              + svi,data = train)
summary(modeloC)
## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + svi, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.69494 -0.46058 -0.04485  0.49149  1.68289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.02278    0.71296  -1.435 0.156362    
## lcavol       0.51999    0.09442   5.507 7.16e-07 ***
## lweight      0.73680    0.20155   3.656 0.000525 ***
## svi          0.53790    0.27093   1.985 0.051458 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7444 on 63 degrees of freedom
## Multiple R-squared:  0.6374, Adjusted R-squared:  0.6202 
## F-statistic: 36.92 on 3 and 63 DF,  p-value: 6.81e-14
predictions <- predict(modeloC,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)
##     RMSE Rsquared 
##    0.633    0.627
axisRange = extendrange(c(test$lpsa,predictions))
plot(test$lpsa,predictions)
abline(0,1,col="blue",lty=2,lwd=2)

Há evidência de relação entre os preditores e a variável alvo?

Sim,como foi possível observar nos gráficos acima, o hiperlano se ajusta muito bem aos dados de teste.

Havendo relação, quão forte é essa relação?

A relação é consideravalmente alta devido a baixa taxa de erro residual. Além disso, é possível analisar visualmente esta forte relação

Que variável parece contribuir mais?

A variável que mais parece contribuir com a predição é a lcavol. Que por sinal também apresentou melhor correlação com a varíavel destino

A relação sugere um modelo de regressão linear?

Sim. O hiperlano se ajusta consideravelmente bem aos dados. No entanto, temos dois outliers que não podem ser explicados pelo modelo devido a quantidade de dados a qual é muito pequena

Abaixo podemos ver o sumário dos 3 modelos utilizando todo o dataset (treino+teste), sendo o último o escolhido

##     RMSE Rsquared 
##    0.682    0.648

##     RMSE Rsquared 
##    0.681    0.648

##     RMSE Rsquared 
##    0.696    0.634