Relatorio_1
Rodolfo Viana
27-03-2016
Câncer de Prostata
Antígeno Prostático Específico (PSA) é uma substância produzida pelas células da glândula prostática. O PSA é encontrado principalmente no sêmen, mas uma pequena quantidade é também encontrada no sangue. A maioria dos homens saudáveis têm níveis menores de 4ng/ml de sangue. Nosso objetivo final será predizer a variável psa a partir dos preditores disponíveis.
Utilizamos nesse problema os dados referentes a exames em pacientes homens, afim de identificar sintomas de câncer de prostata. O dataset é composto por dez variáveis:
Como campos do dataset, temos:
- 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).
library(dplyr)
library(ggplot2)
library(caret)Carregando os dados e dividindo em treino e teste
prostate <- read.delim("/home/rodolfo/Projetos/DataAnalysis/LinearRegression/dados/prostate.data")
prostate.train <- filter(prostate, train == TRUE)
prostate.teste <- filter(prostate, train == FALSE)A princípio iremos fazer uma análise descritiva dos dados e uma análise de correlação das variáveis para identificar possíveis redundâncias.
library(GGally)
ggpairs(select(prostate.train, lcavol, lweight, age, lbph, svi, lcp, gleason, pgg45, lpsa))
O gráfico acima possui diversas informações, para essa parte do problemas focamos apenas na busca de alguma variável que tenha alguma impacto na variável lpsa. Estamos buscando a variável que mais tem correlação com a variável lpsa. Por essa razão escolhemos a variável lcavol por ser a varável que possui maior correlação 0.733.
Analisando o lcavol e o lpsa mais de perto temos
ggplot(prostate.train, aes(lpsa, lcavol)) +
geom_point() +
theme_classic() +
theme(axis.ticks = element_blank(),
legend.position="none")
Podemos notar que essas duas variáveis sugere uma relação linear
Iremos agora construir o nosso primeiro modelo. Esse modelo será bem basico com apenas uma variável e vai servir como base para os outros modelos.
mod_1 <- lm(lpsa ~ lcavol, data = prostate.train)
summary(mod_1)##
## Call:
## lm(formula = lpsa ~ lcavol, data = prostate.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6802 -0.4240 0.1684 0.5392 1.9070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.51630 0.14772 10.264 3.12e-15 ***
## lcavol 0.71264 0.08199 8.692 1.73e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8277 on 65 degrees of freedom
## Multiple R-squared: 0.5375, Adjusted R-squared: 0.5304
## F-statistic: 75.55 on 1 and 65 DF, p-value: 1.733e-12É possível notar que temos alta a significância de que é pouco provável o fato de que não exista nenhuma relação entre as variáveis lpsa e lcavol. Há evidência de que a relação entre essas duas variáveis seja forte, dado que o R-squared igual à 0,53. O que siginifica que a variável escolhidas explica em 53% a variável resposta (lpsa).
Após criado o modelo, iremos prever utilizando os dados de teste. Para entender melhor a diferença do que foi previsto e do valor esperado podemos observar o gráfico abaixo:
predicoes = predict.lm(mod_1, prostate.teste)
toPlot <- data.frame(esperado = prostate.teste$lpsa,
previsto = predicoes)
residuos <- toPlot$esperado - toPlot$previsto
ggplot(toPlot, aes(esperado, previsto)) +
geom_point() +
theme_classic() +
theme(axis.ticks = element_blank(),
legend.position="none")
Em um modelo ideal todos os pontos estariam formando uma linha na diagonal do gráfico. Podemos notar que, considerando que estamos usando apenas uma variável, temos um modelo considerado razoável.
RMSE(toPlot$previsto, toPlot$esperado)## [1] 0.6926317RMSE foi de 0.69
É importante também verificar se os resíduos seguem uma distribuição normal com média 0:
qqnorm(residuos)
qqline(residuos, col = 2,lwd=2,lty=2)
Na maior parte os resíduos seguem uma distribuição normal.
Iremos agora adicionar mais uma variáveis no nosso modelo e verificar se o nosso modelo fica melhor ou pior. Escolhemos as variáveis lcp e svi utilizando o mesmo critério anterior.
mod_3 <- lm(lpsa ~ lcavol + svi, data = prostate.train)
predicoes = predict.lm(mod_3, prostate.teste)
toPlot <- data.frame(esperado = prostate.teste$lpsa,
previsto = predicoes)
residuos <- toPlot$esperado - toPlot$previsto
linear_rmse <- RMSE(toPlot$previsto, toPlot$esperado)
linear_rmse## [1] 0.6213851Com esse novo modelo temos um RMSE mais baixo do que o anterior. O que significa que o nosso modelo errou menos do que o modelo anterior.
Adicionamos mais uma variável ao nosso modelo. Escolhemos as variáveis lcp utilizando o mesmo critério anterior.
mod_2 <- lm(lpsa ~ lcavol + lcp + svi, data = prostate.train)
predicoes = predict.lm(mod_2, prostate.teste)
toPlot <- data.frame(esperado = prostate.teste$lpsa,
previsto = predicoes)
residuos <- toPlot$esperado - toPlot$previsto
RMSE(toPlot$previsto, toPlot$esperado)## [1] 0.6654602Com esse novo modelo temos um RMSE mais alto do que o anterior. O que significa que o nosso modelo errou mais que o anterior. Isso pode ter resultado de uma variável pouco importante para o modelo.
summary(mod_2)##
## Call:
## lm(formula = lpsa ~ lcavol + lcp + svi, data = prostate.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7015 -0.5090 0.1243 0.5160 1.6985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3722 0.1952 7.030 1.77e-09 ***
## lcavol 0.6754 0.1144 5.903 1.55e-07 ***
## lcp -0.1395 0.1103 -1.265 0.2104
## svi 0.7290 0.3296 2.212 0.0306 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8093 on 63 degrees of freedom
## Multiple R-squared: 0.5714, Adjusted R-squared: 0.551
## F-statistic: 28 on 3 and 63 DF, p-value: 1.256e-11Verificamos que a variável lcp tem pouca significancia para o modelo. Por isso, dos três modelos testados o melhor modelo para esse conjunto de dados é o segundo modelo, utilizando como critério o valor do RMSE e o R-squared.
Iremos agora criar um modelo otimizado utilizando Lasso. O Lasso é uma técnica que, além de controlar o overfitting, aplica a seleção de variáveis que melhor explicam a variável resposta. Para criar esse modelo iremos utilizar a biblioteca caret.
Primeiramente vamos criar um modelo no caret com todas as variáveis possíveis.
treino_labels = prostate.train[, 10] ## Classes das instâncias de treino
treino = prostate.train[-c(1,10,11)] ## Exclui variável alvo
lmFit <- train(treino, treino_labels, method= "lm", metric="RMSE")
summary(lmFit)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## 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 *
## lcp -0.206324 0.110516 -1.867 0.06697 .
## gleason -0.029503 0.201136 -0.147 0.88389
## pgg45 0.009465 0.005447 1.738 0.08755 .
## ---
## 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-12Iremos excluir as variáveis com pouca significância (age, lcp, gleason, pgg45). Além disso, vamos realizar o treino utilizando a cross validation.
treino_labels = prostate.train[, 10] ## Classes das instâncias de treino
treino = prostate.train[c(2,3,5,6)] ## Exclui variável alvo
con <- trainControl(method="cv", number=10)
lmFit <- train(treino, treino_labels, method= "lm", metric="RMSE", trControl = con)
test <- prostate.teste[c(2,3,5,6)]
lmfit_prediction <- predict(lmFit, test)
toPlot <- data.frame(esperado = prostate.teste$lpsa,
previsto = lmfit_prediction)
lmfit_rmse <- RMSE(toPlot$previsto, toPlot$esperado)
summary(lmFit)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8709 -0.3903 -0.0172 0.5676 1.4227
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.32592 0.77998 -0.418 0.6775
## lcavol 0.50552 0.09256 5.461 8.85e-07 ***
## lweight 0.53883 0.22071 2.441 0.0175 *
## lbph 0.14001 0.07041 1.988 0.0512 .
## svi 0.67185 0.27323 2.459 0.0167 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7275 on 62 degrees of freedom
## Multiple R-squared: 0.6592, Adjusted R-squared: 0.6372
## F-statistic: 29.98 on 4 and 62 DF, p-value: 6.911e-14Para otimizar ainda mais o nosso modelo iremos utilizar Lasso e a validação cruzada. Temos então o seguinte resultado de treinamento para a validação cruzada
treino_labels = prostate.train[, 10] ## Classes das instâncias de treino
treino = prostate.train[c(2,3,5,6)] ## Exclui variável alvo
con <- trainControl(method="cv", number=10)
lassoFit <- train(treino, treino_labels, method= "lasso", metric="RMSE", trControl = con)
plot(lassoFit)
A função train tentou 10 valores para fraction e achou valor ótimo (RMSE mais baixo).
lassoFit ## The lasso
##
## 67 samples
## 4 predictors
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
##
## Summary of sample sizes: 61, 60, 60, 61, 59, 62, ...
##
## Resampling results across tuning parameters:
##
## fraction RMSE Rsquared RMSE SD Rsquared SD
## 0.1 1.0664612 0.6174823 0.2869788 0.1639269
## 0.5 0.8190214 0.6645981 0.2216896 0.1829279
## 0.9 0.7523764 0.7122103 0.1500642 0.2169375
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 0.9.Depois de criado o modelo iremos agora gerar a previsão
train <- prostate.teste[c(2,3,5,6)]
lasso_prediction <- predict(lassoFit, train)
toPlot <- data.frame(esperado = prostate.teste$lpsa,
previsto = lasso_prediction)
residuos <- toPlot$esperado - toPlot$previsto
lasso_rmse <- RMSE(toPlot$previsto, toPlot$esperado)
lasso_rmse## [1] 0.6662196Comparando o RMSE do melhor modelo utilizando regressão linear x regressão linear com cross validation x lasso
toPlot <- data.frame(RMSE = c(linear_rmse, lmfit_rmse, lasso_rmse),
Modelo = c("Linear Regression", "Linear Regression CV", "Lasso"))
ggplot(toPlot, aes(x=reorder(Modelo, -RMSE), y=RMSE)) +
geom_bar(stat="identity") +
labs(x='Modelo', y='RMSE') +
theme(axis.text.x = element_text(angle = 45, hjust = 1), panel.background=element_blank()) +
coord_flip()
O melhor modelo será aquele que possuir o RMSE mais baixo. Utilizando essa métrica o melhor modelo gerado foi o utilizando linear regression sem o cross validation.