Análise de Dados 2 - UFCG
Diego Pedro
library(corrplot)
library(reshape2)
library(ggplot2)
library(caret)## Loading required package: latticelibrary(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, unionOs 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.5829No 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 variablesggplot(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-12predictions <- 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.505modeloB <- 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-13predictions <- 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.561Defininindo 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-14Analisando 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.627Na 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 ' ' 1Conclusã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(p = leaps.prostate$size,Cp = leaps.prostate$Cp)
ggplot(lm_prediction, aes(x=p,y=Cp)) +
geom_point(shape=1) +
geom_smooth()
No gráfico acima temos no eixo y os coeficientes de Mallows e no eixo x o tamanho do subset de preditores
O coefficiente de Mallows é utilizado para escolher o tamanho ideal de subset de variáveis preditoras. O ideal é quando Cp(coeficiente de Mallow) é próximo ou menor que p(quantidade de parâmetros do subset). No entanto, é preciso ter cuidado para não escolher o menor Cp em relação a p para evitar overffiting. Foi escolhido o subset de tamanho 4 pois é ponto em que o função começa a se estabilizar.
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.
Abaixo o cálculo do BIC
BIC: \(ln(RSS)n - ln(n)n + ln(n)p\)
p é o número de preditores
n é o tamanho da amostra
Ou seja, a medida o tamanho da amostra aumenta BIC tende a decair. No entanto, a medida que o RSS aumenta o BIC tende a subir.
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.2Como 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, 59, 62, 59, 61, 59, ...
## Resampling results across tuning parameters:
##
## fraction RMSE Rsquared RMSE SD Rsquared SD
## 0.1000000 1.0154122 0.5533656 0.3007352 0.2522665
## 0.1888889 0.9145588 0.5525074 0.2673515 0.2470324
## 0.2777778 0.8430578 0.5846421 0.2383203 0.2457313
## 0.3666667 0.7884250 0.6028473 0.2159043 0.2365889
## 0.4555556 0.7530468 0.6154958 0.2034926 0.2070922
## 0.5444444 0.7310697 0.6281613 0.1974430 0.1826297
## 0.6333333 0.7367351 0.6173018 0.1977065 0.1833173
## 0.7222222 0.7378233 0.6131861 0.2034592 0.1850805
## 0.8111111 0.7329498 0.6152073 0.2058487 0.1906779
## 0.9000000 0.7302225 0.6150989 0.2110972 0.2012273
##
## 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.527Aplicando Ridge
Uma outra forma de tratar o overfitting é usando a técnica Ridge.No entanto, esta técnica não faz seleção de preditores
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. O melhor estado será quando lambda tenderá ao infinito.
library(genridge)## Loading required package: carridge.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)Abaixo estão os gráficos conhecidos como Traceplote utilizado para decidir o valor ideal de 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. Claramente podemos observar que as variáveis lcavol e svi são as de maior importância.
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, 60, 60, 59, 60, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared RMSE SD Rsquared SD
## 0.0000000000 0.7324616 0.6691926 0.2366868 0.2953197
## 0.0001000000 0.7324386 0.6692375 0.2367016 0.2953202
## 0.0002371374 0.7324072 0.6692989 0.2367220 0.2953209
## 0.0005623413 0.7323333 0.6694441 0.2367703 0.2953224
## 0.0013335214 0.7321610 0.6697850 0.2368853 0.2953260
## 0.0031622777 0.7317687 0.6705760 0.2371601 0.2953332
## 0.0074989421 0.7309248 0.6723568 0.2378202 0.2953421
## 0.0177827941 0.7293486 0.6760916 0.2394042 0.2953021
## 0.0421696503 0.7273883 0.6827268 0.2430550 0.2947592
## 0.1000000000 0.7282234 0.6906972 0.2505208 0.2910650
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 0.04216965.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.708 0.524Utilizando Regressao Polinomial
Polinômio de Grau 2
modeloP2 <- lm(lpsa ~ poly(lcavol, 2, raw=TRUE)
+poly(lweight, 2, raw=TRUE)
+ poly(age, 2, raw=TRUE)
+ poly(lbph, 2, raw=TRUE)
+ poly(pgg45, 2, raw=TRUE)
+ poly(lcp, 2, raw=TRUE)
+ poly(gleason, 2, raw=TRUE)
+ svi,data = train)
summary(modeloP2)##
## Call:
## lm(formula = lpsa ~ poly(lcavol, 2, raw = TRUE) + poly(lweight,
## 2, raw = TRUE) + poly(age, 2, raw = TRUE) + poly(lbph, 2,
## raw = TRUE) + poly(pgg45, 2, raw = TRUE) + poly(lcp, 2, raw = TRUE) +
## poly(gleason, 2, raw = TRUE) + svi, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.79739 -0.27192 -0.02137 0.40111 1.34662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.384e+00 1.506e+01 0.358 0.72212
## poly(lcavol, 2, raw = TRUE)1 5.465e-01 1.712e-01 3.193 0.00241 **
## poly(lcavol, 2, raw = TRUE)2 2.542e-02 6.882e-02 0.369 0.71341
## poly(lweight, 2, raw = TRUE)1 -1.426e+00 2.401e+00 -0.594 0.55514
## poly(lweight, 2, raw = TRUE)2 3.050e-01 3.287e-01 0.928 0.35772
## poly(age, 2, raw = TRUE)1 -8.036e-02 1.352e-01 -0.594 0.55490
## poly(age, 2, raw = TRUE)2 4.813e-04 1.096e-03 0.439 0.66242
## poly(lbph, 2, raw = TRUE)1 1.718e-01 8.215e-02 2.091 0.04152 *
## poly(lbph, 2, raw = TRUE)2 -1.521e-01 1.018e-01 -1.494 0.14141
## poly(pgg45, 2, raw = TRUE)1 1.708e-02 2.229e-02 0.766 0.44704
## poly(pgg45, 2, raw = TRUE)2 -9.878e-05 2.521e-04 -0.392 0.69681
## poly(lcp, 2, raw = TRUE)1 -1.902e-01 1.190e-01 -1.599 0.11601
## poly(lcp, 2, raw = TRUE)2 -1.045e-01 7.790e-02 -1.342 0.18564
## poly(gleason, 2, raw = TRUE)1 1.758e-01 3.397e+00 0.052 0.95893
## poly(gleason, 2, raw = TRUE)2 -1.403e-02 2.299e-01 -0.061 0.95156
## svi 6.971e-01 3.273e-01 2.130 0.03805 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.714 on 51 degrees of freedom
## Multiple R-squared: 0.7299, Adjusted R-squared: 0.6505
## F-statistic: 9.19 on 15 and 51 DF, p-value: 8.054e-10Podemos observar que o resultado do RMSE foi inferior ao do primeiro modelo usando todas as variáveis a qual foi 0.722
predictions <- predict(modeloP2,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)## RMSE Rsquared
## 0.744 0.478Polinômio de Grau 3
modeloP3 <- lm(lpsa ~ poly(lcavol, 3, raw=TRUE)
+poly(lweight, 3, raw=TRUE)
+ poly(age, 3, raw=TRUE)
+ poly(lbph, 3, raw=TRUE)
+ poly(pgg45, 3, raw=TRUE)
+ poly(lcp, 3, raw=TRUE)
+ poly(gleason, 3, raw=TRUE)
+ svi,data = train)
summary(modeloP3)##
## Call:
## lm(formula = lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight,
## 3, raw = TRUE) + poly(age, 3, raw = TRUE) + poly(lbph, 3,
## raw = TRUE) + poly(pgg45, 3, raw = TRUE) + poly(lcp, 3, raw = TRUE) +
## poly(gleason, 3, raw = TRUE) + svi, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82824 -0.32405 -0.01544 0.30537 1.26523
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.544e+01 2.135e+02 0.400 0.69093
## poly(lcavol, 3, raw = TRUE)1 5.145e-01 1.807e-01 2.848 0.00667 **
## poly(lcavol, 3, raw = TRUE)2 -7.084e-02 1.743e-01 -0.406 0.68647
## poly(lcavol, 3, raw = TRUE)3 3.297e-02 4.912e-02 0.671 0.50553
## poly(lweight, 3, raw = TRUE)1 -1.288e+01 1.842e+01 -0.699 0.48799
## poly(lweight, 3, raw = TRUE)2 3.545e+00 5.174e+00 0.685 0.49679
## poly(lweight, 3, raw = TRUE)3 -3.015e-01 4.775e-01 -0.632 0.53096
## poly(age, 3, raw = TRUE)1 -6.286e-01 1.310e+00 -0.480 0.63368
## poly(age, 3, raw = TRUE)2 9.761e-03 2.182e-02 0.447 0.65687
## poly(age, 3, raw = TRUE)3 -5.089e-05 1.195e-04 -0.426 0.67240
## poly(lbph, 3, raw = TRUE)1 2.847e-01 2.837e-01 1.003 0.32118
## poly(lbph, 3, raw = TRUE)2 -1.076e-01 1.854e-01 -0.580 0.56460
## poly(lbph, 3, raw = TRUE)3 -3.615e-02 1.243e-01 -0.291 0.77249
## poly(pgg45, 3, raw = TRUE)1 -2.721e-02 5.461e-02 -0.498 0.62087
## poly(pgg45, 3, raw = TRUE)2 1.038e-03 1.296e-03 0.801 0.42753
## poly(pgg45, 3, raw = TRUE)3 -8.062e-06 9.364e-06 -0.861 0.39395
## poly(lcp, 3, raw = TRUE)1 -1.979e-01 2.677e-01 -0.739 0.46365
## poly(lcp, 3, raw = TRUE)2 -1.168e-01 1.476e-01 -0.791 0.43297
## poly(lcp, 3, raw = TRUE)3 3.343e-03 8.167e-02 0.041 0.96754
## poly(gleason, 3, raw = TRUE)1 -2.421e+01 9.003e+01 -0.269 0.78928
## poly(gleason, 3, raw = TRUE)2 3.442e+00 1.238e+01 0.278 0.78234
## poly(gleason, 3, raw = TRUE)3 -1.591e-01 5.612e-01 -0.284 0.77808
## svi 7.651e-01 3.638e-01 2.103 0.04118 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7459 on 44 degrees of freedom
## Multiple R-squared: 0.7457, Adjusted R-squared: 0.6186
## F-statistic: 5.866 on 22 and 44 DF, p-value: 3.308e-07Podemos observar que Polinômio de grau 3 mostrou piora
predictions <- predict(modeloP3,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)## RMSE Rsquared
## 0.751 0.469Polinômio de Grau 3 com seleção de preditores (realizado anteriormente)
modeloP3B <- lm(lpsa ~ poly(lcavol, 3, raw=TRUE)
+poly(lweight, 3, raw=TRUE)
+ svi,data = train)
summary(modeloP3B)##
## Call:
## lm(formula = lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight,
## 3, raw = TRUE) + svi, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.73161 -0.49381 0.04038 0.55115 1.69921
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.17870 18.66262 0.760 0.4504
## poly(lcavol, 3, raw = TRUE)1 0.65670 0.14849 4.423 4.26e-05 ***
## poly(lcavol, 3, raw = TRUE)2 -0.20157 0.15445 -1.305 0.1969
## poly(lcavol, 3, raw = TRUE)3 0.03944 0.04111 0.959 0.3413
## poly(lweight, 3, raw = TRUE)1 -12.34230 15.87393 -0.778 0.4400
## poly(lweight, 3, raw = TRUE)2 3.71209 4.44096 0.836 0.4066
## poly(lweight, 3, raw = TRUE)3 -0.34468 0.40869 -0.843 0.4024
## svi 0.67284 0.30162 2.231 0.0295 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7495 on 59 degrees of freedom
## Multiple R-squared: 0.6557, Adjusted R-squared: 0.6149
## F-statistic: 16.06 on 7 and 59 DF, p-value: 1.299e-11Já com a seleção de variáveis obtemos uma melhora considerável no RMSE. No entanto, ainda é inferior ao do modelo de polinômio 1.
predictions <- predict(modeloP3B,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)## RMSE Rsquared
## 0.646 0.613Por fim utilizando o stepwise automaticamente para o modelo de Polinômio 3
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectstep <- stepAIC(modeloP3, direction="both")## Start: AIC=-21.46
## lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight, 3, raw = TRUE) +
## poly(age, 3, raw = TRUE) + poly(lbph, 3, raw = TRUE) + poly(pgg45,
## 3, raw = TRUE) + poly(lcp, 3, raw = TRUE) + poly(gleason,
## 3, raw = TRUE) + svi
##
## Df Sum of Sq RSS AIC
## - poly(gleason, 3, raw = TRUE) 3 0.4157 24.896 -26.3288
## - poly(age, 3, raw = TRUE) 3 0.6795 25.160 -25.6228
## - poly(pgg45, 3, raw = TRUE) 3 1.6156 26.096 -23.1752
## <none> 24.480 -21.4571
## - poly(lcp, 3, raw = TRUE) 3 2.4245 26.905 -21.1299
## - poly(lbph, 3, raw = TRUE) 3 3.2031 27.683 -19.2185
## - poly(lweight, 3, raw = TRUE) 3 3.9024 28.383 -17.5472
## - svi 1 2.4616 26.942 -17.0376
## - poly(lcavol, 3, raw = TRUE) 3 12.6210 37.101 0.4003
##
## Step: AIC=-26.33
## lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight, 3, raw = TRUE) +
## poly(age, 3, raw = TRUE) + poly(lbph, 3, raw = TRUE) + poly(pgg45,
## 3, raw = TRUE) + poly(lcp, 3, raw = TRUE) + svi
##
## Df Sum of Sq RSS AIC
## - poly(age, 3, raw = TRUE) 3 0.6864 25.582 -30.5066
## <none> 24.896 -26.3288
## - poly(pgg45, 3, raw = TRUE) 3 2.4195 27.316 -26.1148
## - poly(lcp, 3, raw = TRUE) 3 2.6636 27.560 -25.5187
## - poly(lbph, 3, raw = TRUE) 3 3.2050 28.101 -24.2153
## - poly(lweight, 3, raw = TRUE) 3 3.8478 28.744 -22.7001
## - svi 1 2.3560 27.252 -22.2706
## + poly(gleason, 3, raw = TRUE) 3 0.4157 24.480 -21.4571
## - poly(lcavol, 3, raw = TRUE) 3 14.1685 39.065 -2.1449
##
## Step: AIC=-30.51
## lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight, 3, raw = TRUE) +
## poly(lbph, 3, raw = TRUE) + poly(pgg45, 3, raw = TRUE) +
## poly(lcp, 3, raw = TRUE) + svi
##
## Df Sum of Sq RSS AIC
## - poly(pgg45, 3, raw = TRUE) 3 2.3397 27.922 -30.6431
## <none> 25.582 -30.5066
## - poly(lcp, 3, raw = TRUE) 3 2.5845 28.167 -30.0583
## - poly(lbph, 3, raw = TRUE) 3 3.1578 28.740 -28.7083
## - poly(lweight, 3, raw = TRUE) 3 3.6394 29.222 -27.5950
## - svi 1 2.2572 27.840 -26.8414
## + poly(age, 3, raw = TRUE) 3 0.6864 24.896 -26.3288
## + poly(gleason, 3, raw = TRUE) 3 0.4227 25.160 -25.6228
## - poly(lcavol, 3, raw = TRUE) 3 14.4645 40.047 -6.4808
##
## Step: AIC=-30.64
## lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight, 3, raw = TRUE) +
## poly(lbph, 3, raw = TRUE) + poly(lcp, 3, raw = TRUE) + svi
##
## Df Sum of Sq RSS AIC
## - poly(lcp, 3, raw = TRUE) 3 1.9839 29.906 -32.044
## <none> 27.922 -30.643
## + poly(pgg45, 3, raw = TRUE) 3 2.3397 25.582 -30.507
## - poly(lbph, 3, raw = TRUE) 3 3.5752 31.497 -28.571
## + poly(gleason, 3, raw = TRUE) 3 0.9796 26.943 -27.036
## - poly(lweight, 3, raw = TRUE) 3 4.5418 32.464 -26.545
## + poly(age, 3, raw = TRUE) 3 0.6066 27.316 -26.115
## - svi 1 2.9909 30.913 -25.825
## - poly(lcavol, 3, raw = TRUE) 3 15.5624 43.485 -6.963
##
## Step: AIC=-32.04
## lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight, 3, raw = TRUE) +
## poly(lbph, 3, raw = TRUE) + svi
##
## Df Sum of Sq RSS AIC
## <none> 29.906 -32.044
## - poly(lbph, 3, raw = TRUE) 3 3.2389 33.145 -31.155
## + poly(lcp, 3, raw = TRUE) 3 1.9839 27.922 -30.643
## + poly(pgg45, 3, raw = TRUE) 3 1.7391 28.167 -30.058
## - svi 1 2.2146 32.121 -29.258
## - poly(lweight, 3, raw = TRUE) 3 4.4465 34.353 -28.757
## + poly(gleason, 3, raw = TRUE) 3 1.1817 28.724 -28.745
## + poly(age, 3, raw = TRUE) 3 0.8488 29.057 -27.973
## - poly(lcavol, 3, raw = TRUE) 3 16.5148 46.421 -8.585summary(step)##
## Call:
## lm(formula = lpsa ~ poly(lcavol, 3, raw = TRUE) + poly(lweight,
## 3, raw = TRUE) + poly(lbph, 3, raw = TRUE) + svi, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.01884 -0.34826 -0.01174 0.44216 1.45512
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.69761 18.61641 0.736 0.464934
## poly(lcavol, 3, raw = TRUE)1 0.58148 0.15034 3.868 0.000289 ***
## poly(lcavol, 3, raw = TRUE)2 -0.19410 0.15142 -1.282 0.205165
## poly(lcavol, 3, raw = TRUE)3 0.04942 0.04089 1.209 0.231877
## poly(lweight, 3, raw = TRUE)1 -10.76052 15.83994 -0.679 0.499728
## poly(lweight, 3, raw = TRUE)2 3.01593 4.43849 0.679 0.499623
## poly(lweight, 3, raw = TRUE)3 -0.25865 0.40882 -0.633 0.529520
## poly(lbph, 3, raw = TRUE)1 0.27306 0.24375 1.120 0.267406
## poly(lbph, 3, raw = TRUE)2 -0.10615 0.15531 -0.683 0.497125
## poly(lbph, 3, raw = TRUE)3 -0.04134 0.10343 -0.400 0.690890
## svi 0.62359 0.30622 2.036 0.046449 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7308 on 56 degrees of freedom
## Multiple R-squared: 0.6894, Adjusted R-squared: 0.6339
## F-statistic: 12.43 on 10 and 56 DF, p-value: 5.237e-11Podemos observar que o stepwise incluiu a variável lbph no modelo, diferentemente da seleção manual.
predictions <- predict(step,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)## RMSE Rsquared
## 0.677 0.563Por fim utilizando o stepwise automaticamente para o modelo de Regressão Linear
library(MASS)
step <- stepAIC(modeloA, direction="both")## Start: AIC=-37.13
## lpsa ~ lcavol + lweight + age + lbph + svi + pgg45 + lcp + gleason
##
## Df Sum of Sq RSS AIC
## - gleason 1 0.0109 29.437 -39.103
## <none> 29.426 -37.128
## - age 1 0.9886 30.415 -36.914
## - pgg45 1 1.5322 30.959 -35.727
## - lcp 1 1.7683 31.195 -35.218
## - lbph 1 2.1443 31.571 -34.415
## - svi 1 3.0934 32.520 -32.430
## - lweight 1 3.8390 33.265 -30.912
## - lcavol 1 14.6102 44.037 -12.118
##
## Step: AIC=-39.1
## lpsa ~ lcavol + lweight + age + lbph + svi + pgg45 + lcp
##
## Df Sum of Sq RSS AIC
## <none> 29.437 -39.103
## - age 1 1.1025 30.540 -38.639
## - lcp 1 1.7583 31.196 -37.216
## + gleason 1 0.0109 29.426 -37.128
## - lbph 1 2.1354 31.573 -36.411
## - pgg45 1 2.3755 31.813 -35.903
## - svi 1 3.1665 32.604 -34.258
## - lweight 1 4.0048 33.442 -32.557
## - lcavol 1 14.8873 44.325 -13.681summary(step)##
## Call:
## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + pgg45 +
## lcp, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65425 -0.34471 -0.05615 0.44380 1.48952
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.259062 1.025170 0.253 0.8014
## lcavol 0.573930 0.105069 5.462 9.88e-07 ***
## lweight 0.619209 0.218560 2.833 0.0063 **
## age -0.019480 0.013105 -1.486 0.1425
## lbph 0.144426 0.069812 2.069 0.0430 *
## svi 0.741781 0.294451 2.519 0.0145 *
## pgg45 0.008945 0.004099 2.182 0.0331 *
## lcp -0.205417 0.109424 -1.877 0.0654 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7064 on 59 degrees of freedom
## Multiple R-squared: 0.6943, Adjusted R-squared: 0.658
## F-statistic: 19.14 on 7 and 59 DF, p-value: 4.496e-13No caso do modelo acima o stepwise do R não selecionou nenhuma varíavel
predictions <- predict(step,test)
lm_prediction <- data.frame(pred = predictions, obs = test$lpsa)
round(defaultSummary(lm_prediction), digits = 3)## RMSE Rsquared
## 0.719 0.509Conclusão
Podemos concluir que as técnicas utilizadas acima apresentam desempenho inferior ao modelo do primeiro problema e que o modelo polinomial não apresenta melhoras em relação ao modelo de de regressão linear de polinômio 1
Vale salientar que a complexidade do último modelo apresentado é menor, logo a variância é menor, mas o bias tende a ser maior.