Análise de Dados 2 - UFCG

Diego Pedro

library(corrplot)
library(reshape2)
library(ggplot2)
library(caret)

## Loading required package: lattice

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Os dados concistem 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")

Sumarizando os Dados

summary(prostate.data)

##        X          lcavol           lweight           age       
##  Min.   : 1   Min.   :-1.3471   Min.   :2.375   Min.   :41.00  
##  1st Qu.:25   1st Qu.: 0.5128   1st Qu.:3.376   1st Qu.:60.00  
##  Median :49   Median : 1.4469   Median :3.623   Median :65.00  
##  Mean   :49   Mean   : 1.3500   Mean   :3.629   Mean   :63.87  
##  3rd Qu.:73   3rd Qu.: 2.1270   3rd Qu.:3.876   3rd Qu.:68.00  
##  Max.   :97   Max.   : 3.8210   Max.   :4.780   Max.   :79.00  
##       lbph              svi              lcp             gleason     
##  Min.   :-1.3863   Min.   :0.0000   Min.   :-1.3863   Min.   :6.000  
##  1st Qu.:-1.3863   1st Qu.:0.0000   1st Qu.:-1.3863   1st Qu.:6.000  
##  Median : 0.3001   Median :0.0000   Median :-0.7985   Median :7.000  
##  Mean   : 0.1004   Mean   :0.2165   Mean   :-0.1794   Mean   :6.753  
##  3rd Qu.: 1.5581   3rd Qu.:0.0000   3rd Qu.: 1.1787   3rd Qu.:7.000  
##  Max.   : 2.3263   Max.   :1.0000   Max.   : 2.9042   Max.   :9.000  
##      pgg45             lpsa           train        
##  Min.   :  0.00   Min.   :-0.4308   Mode :logical  
##  1st Qu.:  0.00   1st Qu.: 1.7317   FALSE:30       
##  Median : 15.00   Median : 2.5915   TRUE :67       
##  Mean   : 24.38   Mean   : 2.4784   NA's :0        
##  3rd Qu.: 40.00   3rd Qu.: 3.0564                  
##  Max.   :100.00   Max.   : 5.5829

No gráfico abaixo é apresentado as distribuições das variáveis onde no eixo x representa o valor correspondente da variável e no eixo y a quantidade daquele determinado valor que apareceu nos dados. Min: Menor valor da varíavel

1st Qu. Valor referente ao maior valor dos 25% menores valores

Median Valor que está no centro dos valores, ou seja, a quantidade de items com valores superiores a Median é igual ao número de items com valores inferiores a Median.

Mean Média dos valores

3rd Qu Valor referente ao menor valor dos 25% maiores valores.

Max Valor máximo

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, 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.

Selecionando Features

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

Defininindo o melhor modelo

Por fim, utilizando a técnica Stepwise chegamos ao nosso modelo final

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

Analisando os Resultados

predictions <- predict(modeloC,test)
axisRange = extendrange(c(test$lpsa,predictions))
plot(test$lpsa,predictions)
abline(0,1,col="blue",lty=2,lwd=2)

Uma outra forma de ver este gráfico pode ser visto abaixo

lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
ggplot(lm_prediction, aes(x=obs, y=predictions)) +
    geom_point(shape=1) +    
    geom_smooth()

Podemos observar que existe um “buraco” nos dados devido aos outilers o que acaba acarretando um intervalo de confiança muito alto.

Erro do nosso modelo final, a qual tem o menor RMSE

round(defaultSummary(lm_prediction), digits = 3)

##     RMSE Rsquared 
##    0.633    0.627

Na tabela da Anova abaixo podemos observar que a variável lcavol têm grande importância na explicação da variação dos dados.

anova(modeloC)

## Analysis of Variance Table
## 
## Response: lpsa
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## lcavol     1 51.753  51.753 93.4013 4.664e-14 ***
## lweight    1  7.437   7.437 13.4215 0.0005117 ***
## svi        1  2.184   2.184  3.9418 0.0514583 .  
## Residuals 63 34.908   0.554                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusão

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 e o erro é eplicação da variação dos dados é significativa

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 responsável por exiplicar mais de 51% da variação. Além disso foi a que 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.Caso tivessémos uma quantidade considerável de outilers em relação ao modelo poderíamos considerar que o modelo de regressão linear não seria aplicável,

Parte 2

train <- prostate.data[prostate.data$train==TRUE,]
test <- prostate.data[prostate.data$train==FALSE,]

Selecionado o menor subconjunto de variáveis preditoras: Leaps and Bounds

As vezes quando temos uma quantidade de variáveis preditoras muito grande utilizar os métodos como forward e stepwised pode não ser a melhor alternativa. A técnica leaps and bounds fornece uma estratégia eficiente para encontrar o menor subconjunto de variáveis preditoras que oferecem o menor RSS.

Esta estratégia é bastante útil para enfrentar o problema do trade-off de bias e variância. Se aumentarmos a complexidade de nosso modelo, aumento do número de variáveis, o bias diminui, mas a variância aumenta. Por outra lado, se diminuirmos a complexidade de nosso modelo a variância diminui, mas o bias aumenta.

library(leaps)

predictors <- train[,2:8]
response <- train$lpsa
leaps.prostate <- leaps(predictors,response)

lm_prediction <- data.frame(size = leaps.prostate$size,mallows_cp = leaps.prostate$Cp)
ggplot(lm_prediction, aes(x=size,y=mallows_cp)) +
  geom_point(shape=1) +    
  geom_smooth()

No gráfico acima cada ponto representa um modelo de regressão treinando com um número K de variáveis preditoras e no eixo y o coeficiente Mallows. No problema anterior foi utilizado a técnica stepwise com 3 variáveis. Acima podemos vemos que 4 variáveis seria o ideal.

Utilizando a métrica BIC O código abaixo utiliza a mesma técnica executada acima. No entanto, utiliza uma métrica diferente.

leaps.prostate <- leaps(predictors,response)
subsets <- regsubsets(lpsa ~ lcavol
                      + lweight
                      + age
                      + lbph
                      + svi
                      + pgg45
                      + lcp
                      + gleason,data = train)

plot(subsets)

Cada linha do gráfico acima representa um modelo. Como podemos ver que os melhores modelos, segundo a métrica BIC, usam as variáveis lcavol,lweight,lbph e svi em 3 combinações diferentes. Pode perceber que a medida que a quantidade de variáveis aumenta a métrica BIC começa a decair.

Aplicando Cross Validation

Cross validation é uma técnica utilizada para evitar overfitting do treino, ou seja, enviesamento do algoritmo em relação aos dados de treino. Os dados são divididos em pastas. Uma certa porcentagem de pastas são separadas para teste e a outra parte para treino. Este processo é repetido n vezes, sendo n o número de pastas.

fitControl <- trainControl(method='cv', number = 10)

Aplicando o Lasso

Assim como a técnica Leaps and Bounds, o lasso é um method utilizado para tratar overfitting.

lasso.fit <- train(train$lpsa ~ ., data=select(train,lcavol,lweight,pgg45,age,lbph,svi,lcp,gleason), 
                   method='lasso', 
                   metric="RMSE",
                   tuneLength = 10,
                   trControl=fitControl)

## Loading required package: elasticnet

## Loading required package: lars

## Loaded lars 1.2

Como podemos observar, a medida que eu vou aumentando x o RSME decresce

Plot do Lasso

plot(lasso.fit)

Valores Estimados

lasso.fit

## The lasso 
## 
## 67 samples
##  7 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 60, 61, 60, 61, 62, 60, ... 
## Resampling results across tuning parameters:
## 
##   fraction   RMSE       Rsquared   RMSE SD    Rsquared SD
##   0.1000000  1.0085605  0.5010286  0.2383444  0.3456143  
##   0.1888889  0.9186511  0.5029287  0.1822082  0.3468288  
##   0.2777778  0.8582974  0.5207140  0.1609315  0.3415239  
##   0.3666667  0.8170487  0.5314787  0.1460253  0.3415542  
##   0.4555556  0.7858215  0.5484902  0.1414614  0.3463855  
##   0.5444444  0.7642861  0.5668594  0.1495999  0.3434472  
##   0.6333333  0.7690069  0.5747485  0.1601749  0.3423148  
##   0.7222222  0.7654422  0.5837152  0.1656417  0.3409328  
##   0.8111111  0.7636216  0.5905823  0.1703462  0.3385785  
##   0.9000000  0.7620444  0.5950673  0.1762522  0.3365696  
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was fraction = 0.9.

lasso_prediction <- predict(lasso.fit, select(test,lcavol,lweight,pgg45,age,lbph,svi,lcp,gleason))

lasso_prediction <- data.frame(pred = lasso_prediction, obs = test$lpsa)
round(defaultSummary(lasso_prediction), digits = 3)

##     RMSE Rsquared 
##    0.706    0.527

Aplicando Ridge

Uma outra forma de tratar o overfitting é usando a técnica Ridge.

A técnica Ridge penaliza coeficientes (com execeção do intercept) que têm peso muito alto utilizando a soma dos quadrados dos parâmetros, uma vez que coeficientes com peso muito alto pode prejudicar a generalização.Assim como no Lasso, este algoritmo utiliza um valor lambda a qual é setado automaticamente no código abaixo.

library(genridge)

## Loading required package: car

ridge.fit <- lm(lpsa ~ lcavol
        + lweight
        + age
        + lbph
        + svi
        + pgg45
        + lcp
        + gleason,data = train)
y <- train$lpsa
X0 <- model.matrix(ridge.fit)[,-1]
lambda <- seq(0.001, 100, .01)
aridge0 <- ridge(y, X0, lambda=lambda)
traceplot(aridge0)

traceplot(aridge0, X="df")

Como podemos ver nos dois gráficos acima, conforme a constante de ridge diminui a importância das variáveis aumentam.

Resultados Ridge

Abaixo o gráfico que mostra o desempenho do Ridge utilizando Cross Validation

ridge.fit <- train(lpsa ~ lcavol
                   + lweight
                   + age
                   + lbph
                   + svi
                   + pgg45
                   + lcp
                   + gleason,data = train, 
                   method='ridge', 
                   metric="RMSE",
                   tuneLength = 10,
                   trControl=fitControl)
plot(ridge.fit)

Como podemos observar no gráfico acima. Conforme o peso dos coeficientes diminuem o RSME também. No entanto, existe um ponto que ocorre um aumento repentido do RSME.

Abaixo podemos ver este comportamento em números

ridge.fit

## Ridge Regression 
## 
## 67 samples
## 10 predictors
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 60, 60, 61, 60, 61, 59, ... 
## Resampling results across tuning parameters:
## 
##   lambda        RMSE       Rsquared   RMSE SD    Rsquared SD
##   0.0000000000  0.7357177  0.6928691  0.1338084  0.2686691  
##   0.0001000000  0.7357211  0.6928720  0.1338309  0.2686627  
##   0.0002371374  0.7357259  0.6928759  0.1338617  0.2686538  
##   0.0005623413  0.7357374  0.6928850  0.1339352  0.2686328  
##   0.0013335214  0.7357658  0.6929053  0.1341106  0.2685833  
##   0.0031622777  0.7358398  0.6929472  0.1345333  0.2684675  
##   0.0074989421  0.7360491  0.6930131  0.1355699  0.2682002  
##   0.0177827941  0.7367020  0.6930093  0.1381850  0.2675990  
##   0.0421696503  0.7388165  0.6923669  0.1449237  0.2662837  
##   0.1000000000  0.7450813  0.6892072  0.1616652  0.2634503  
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was lambda = 0.

obs: Uma vez que as técnicas acima obtem as amostras de forma randômica os resultados podem variar.

Como podemos observar o RSME utilizando ridge apresentou um resultado bem próximo do lasso.

ridge_prediction <- predict(ridge.fit,test[,2:9])
ridge_prediction <- data.frame(pred = ridge_prediction, obs = test$lpsa)
round(defaultSummary(ridge_prediction), digits = 3)

##     RMSE Rsquared 
##    0.722    0.505

Utilizando Regressao Polinomial

train$lcavol2 <- train$lcavol^2
train$lcavol3 <- train$lcavol2^2
train$lweight2 <- train$lweight^2
train$lweight3 <- train$lweight2^2

test$lcavol2 <- test$lcavol^2
test$lcavol3 <- test$lcavol2^2
test$lweight2 <- test$lweight^2
test$lweight3 <- test$lweight2^2

modeloC <- lm(lpsa ~ lcavol + lcavol2 + lcavol3
              +lweight + lweight2 + lweight3
              + svi ,data = train)
summary(modeloC)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lcavol2 + lcavol3 + lweight + lweight2 + 
##     lweight3 + svi, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.70323 -0.49863  0.00445  0.57276  1.72961 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.110466  14.345557   0.705   0.4837    
## lcavol       0.716329   0.161914   4.424 4.23e-05 ***
## lcavol2     -0.188873   0.121616  -1.553   0.1258    
## lcavol3      0.009256   0.007930   1.167   0.2479    
## lweight     -7.760657  10.801919  -0.718   0.4753    
## lweight2     1.808495   2.259543   0.800   0.4267    
## lweight3    -0.022964   0.028189  -0.815   0.4185    
## svi          0.691851   0.300989   2.299   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7473 on 59 degrees of freedom
## Multiple R-squared:  0.6578, Adjusted R-squared:  0.6172 
## F-statistic:  16.2 on 7 and 59 DF,  p-value: 1.097e-11

predictions <- predict(modeloC,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)

##     RMSE Rsquared 
##    0.647    0.610

Executando o Stepwise manualmente Apos remover as variaveis lweight

modeloC <- lm(lpsa ~ lcavol + lcavol2 + lcavol3
              + svi ,data = train)
summary(modeloC)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lcavol2 + lcavol3 + svi, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7844 -0.5607  0.1426  0.5762  1.7359 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.654371   0.169105   9.783 3.46e-14 ***
## lcavol       0.810491   0.173129   4.681 1.60e-05 ***
## lcavol2     -0.195837   0.130381  -1.502   0.1382    
## lcavol3      0.009849   0.008494   1.159   0.2507    
## svi          0.698196   0.326254   2.140   0.0363 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8109 on 62 degrees of freedom
## Multiple R-squared:  0.5766, Adjusted R-squared:  0.5492 
## F-statistic: 21.11 on 4 and 62 DF,  p-value: 5.085e-11

predictions <- predict(modeloC,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)

##     RMSE Rsquared 
##    0.630    0.637

Por fim utilizando o stepwise automaticamente

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

step <- stepAIC(modeloC, direction="both")

## Start:  AIC=-23.28
## lpsa ~ lcavol + lcavol2 + lcavol3 + svi
## 
##           Df Sum of Sq    RSS      AIC
## - lcavol3  1    0.8840 41.653 -23.8463
## <none>                 40.769 -23.2835
## - lcavol2  1    1.4835 42.253 -22.8888
## - svi      1    3.0115 43.781 -20.5087
## - lcavol   1   14.4110 55.180  -5.0041
## 
## Step:  AIC=-23.85
## lpsa ~ lcavol + lcavol2 + svi
## 
##           Df Sum of Sq    RSS      AIC
## - lcavol2  1    0.6595 42.313 -24.7937
## <none>                 41.653 -23.8463
## + lcavol3  1    0.8840 40.769 -23.2835
## - svi      1    2.8646 44.518 -21.3900
## - lcavol   1   14.0861 55.739  -6.3287
## 
## Step:  AIC=-24.79
## lpsa ~ lcavol + svi
## 
##           Df Sum of Sq    RSS      AIC
## <none>                 42.313 -24.7937
## + lcavol2  1    0.6595 41.653 -23.8463
## - svi      1    2.2160 44.529 -23.3736
## + lcavol3  1    0.0600 42.253 -22.8888
## - lcavol   1   24.1098 66.422   3.4199

summary(step)

## 
## Call:
## lm(formula = lpsa ~ lcavol + svi, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6200 -0.5176  0.1609  0.5256  1.6946 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.5376     0.1456  10.561 1.17e-15 ***
## lcavol        0.6040     0.1000   6.039 8.70e-08 ***
## svi           0.5418     0.2959   1.831   0.0718 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8131 on 64 degrees of freedom
## Multiple R-squared:  0.5605, Adjusted R-squared:  0.5468 
## F-statistic: 40.82 on 2 and 64 DF,  p-value: 3.747e-12

predictions <- predict(step,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)

##     RMSE Rsquared 
##    0.621    0.641

O resultado acima apresentou um resultado melhor do que a regressão do problema anterior.

Conclusão

Podemos concluir que as técnicas utilizadas acima apresentam desempenho inferior ao modelo normal e que o modelo polinomial apresenta melhores resultados se utilizados 4 variaveis de predicao. No entanto, apresenta um resultado inferior ao utilizado apenas 2 variáveis.

Vale salientar que a complexidade do último modelo apresentado é menor, logo a variância é menor, mas o bias tende a ser maior. Este resultado acima pode ser devido a um overfitting de teste uma vez que a base de dados é pequena. Uma alternativa seria utilizar dados de validação dos quais não fazem parte nem da base de teste como também de treino