gls_conejomalo
Mateo Linares
18/4/2020
## Installing package into 'C:/Users/teoli/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)## package 'jpeg' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\teoli\AppData\Local\Temp\Rtmp4Qu6PB\downloaded_packages
Datos
library(readxl)
conejomalo <- read_excel("C:/Users/teoli/Documents/ECONOMETRIA_AVANZADA_2020/conejomalo.xlsx")
View(conejomalo) Grafico de popularidad
library(ggplot2)
pop <- ggplot(conejomalo,aes(x=track_name, y=Rep))+geom_col()+theme(axis.text.x = element_text(angle = 90))Modelo lineal
bb_lm <- lm(Rep ~ danceability+duration_ms+energy+loudness+mode+speechiness+tempo+valence,data=conejomalo)
summary(bb_lm,correlation = T)##
## Call:
## lm(formula = Rep ~ danceability + duration_ms + energy + loudness +
## mode + speechiness + tempo + valence, data = conejomalo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -159.48 -71.45 -16.76 39.97 562.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.105e+02 4.424e+02 1.606 0.1175
## danceability 1.701e+02 2.228e+02 0.763 0.4505
## duration_ms -2.325e-04 5.741e-04 -0.405 0.6880
## energy -5.348e+02 2.777e+02 -1.926 0.0626 .
## loudness 4.023e+01 2.097e+01 1.919 0.0634 .
## mode -3.145e+01 4.178e+01 -0.753 0.4568
## speechiness 6.865e+02 2.639e+02 2.601 0.0136 *
## tempo -3.679e-01 6.596e-01 -0.558 0.5807
## valence -1.774e+02 1.028e+02 -1.726 0.0934 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 133.8 on 34 degrees of freedom
## Multiple R-squared: 0.2672, Adjusted R-squared: 0.09476
## F-statistic: 1.55 on 8 and 34 DF, p-value: 0.1771
##
## Correlation of Coefficients:
## (Intercept) danceability duration_ms energy loudness mode
## danceability -0.62
## duration_ms -0.55 0.28
## energy -0.78 0.26 0.23
## loudness 0.75 -0.18 -0.21 -0.69
## mode 0.04 -0.12 -0.05 0.01 0.08
## speechiness 0.10 0.23 -0.25 -0.18 0.21 -0.13
## tempo -0.36 0.02 0.10 0.17 -0.22 0.03
## valence -0.16 -0.08 0.14 -0.14 -0.27 -0.02
## speechiness tempo
## danceability
## duration_ms
## energy
## loudness
## mode
## speechiness
## tempo -0.11
## valence -0.21 0.18pop
Segundo modelo
bb_lm2 <- lm(Rep ~ energy+loudness+speechiness+valence,data=conejomalo)
summary(bb_lm2,correlation = T) ##
## Call:
## lm(formula = Rep ~ energy + loudness + speechiness + valence,
## data = conejomalo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -145.75 -80.92 -21.48 32.10 600.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 726.44 258.56 2.810 0.0078 **
## energy -533.61 252.73 -2.111 0.0414 *
## loudness 39.63 19.24 2.059 0.0464 *
## speechiness 553.66 231.59 2.391 0.0219 *
## valence -153.15 96.41 -1.588 0.1205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 129.7 on 38 degrees of freedom
## Multiple R-squared: 0.2305, Adjusted R-squared: 0.1495
## F-statistic: 2.846 on 4 and 38 DF, p-value: 0.03704
##
## Correlation of Coefficients:
## (Intercept) energy loudness speechiness
## energy -0.89
## loudness 0.90 -0.66
## speechiness 0.14 -0.19 0.22
## valence -0.16 -0.19 -0.25 -0.14Modelo GLS con cholesky
cor(bb_lm2$res[-1],bb_lm2$res[-43])## [1] 0.2614102x2 <- model.matrix(bb_lm2)
Sigma2<- diag(43)
Sigma2<- 0.2614102^abs(row(Sigma2)-col(Sigma2))
Sigi2<- solve(Sigma2)
xtxi2<- solve(t(x2)%*%Sigi2%*%x2)
beta2<- xtxi2%*% t(x2)%*% Sigi2%*% conejomalo$Rep
beta2## [,1]
## (Intercept) 604.93349
## energy -409.50397
## loudness 31.24595
## speechiness 530.44943
## valence -159.48547sm<- chol(Sigma2)
smi<-solve(t(sm))
sx<-smi%*%x2
sy<-smi%*%conejomalo$Rep
lm(sy~sx-1)$coef## sx(Intercept) sxenergy sxloudness sxspeechiness sxvalence
## 604.93349 -409.50397 31.24595 530.44943 -159.48547Chol_BB<-lm(sy~sx-1)
cor(Chol_BB$res[-1],Chol_BB$res[-43])## [1] 0.2145829GLS con nlme
library(nlme)## Warning: package 'nlme' was built under R version 3.6.3GLS_BB<-gls(Rep ~ energy+loudness+speechiness+valence,correlation =corLin(form=~tempo),data=conejomalo)
summary(GLS_BB)## Generalized least squares fit by REML
## Model: Rep ~ energy + loudness + speechiness + valence
## Data: conejomalo
## AIC BIC logLik
## 496.9548 508.4179 -241.4774
##
## Correlation Structure: Linear spatial correlation
## Formula: ~tempo
## Parameter estimate(s):
## range
## 0.04404024
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 786.2087 258.19385 3.045033 0.0042
## energy -618.9415 251.19951 -2.463944 0.0184
## loudness 42.4168 19.18163 2.211321 0.0331
## speechiness 524.1234 215.80988 2.428635 0.0200
## valence -114.6742 81.55199 -1.406148 0.1678
##
## Correlation:
## (Intr) energy lodnss spchns
## energy -0.899
## loudness 0.894 -0.668
## speechiness 0.246 -0.346 0.280
## valence -0.240 -0.084 -0.301 0.013
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.1666893 -0.6208779 -0.1983574 0.1862973 4.6402076
##
## Residual standard error: 129.3403
## Degrees of freedom: 43 total; 38 residualintervals(GLS_BB)## Approximate 95% confidence intervals
##
## Coefficients:
## lower est. upper
## (Intercept) 263.52260 786.20873 1308.89486
## energy -1127.46831 -618.94149 -110.41467
## loudness 3.58557 42.41676 81.24795
## speechiness 87.23916 524.12341 961.00767
## valence -279.76753 -114.67416 50.41921
## attr(,"label")
## [1] "Coefficients:"
##
## Correlation structure:
## lower est. upper
## range 0.0280745 0.04404024 0.07691563
## attr(,"label")
## [1] "Correlation structure:"
##
## Residual standard error:
## lower est. upper
## 103.2353 129.3403 162.0464GLS con errores robustos y con Sandwich
library(robustbase)## Warning: package 'robustbase' was built under R version 3.6.3library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.6.3## -- Attaching packages ---------------- tidyverse 1.3.0 --## v tibble 2.1.3 v dplyr 0.8.5
## v tidyr 1.0.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## v purrr 0.3.3## Warning: package 'tidyr' was built under R version 3.6.3## Warning: package 'readr' was built under R version 3.6.3## Warning: package 'purrr' was built under R version 3.6.3## Warning: package 'dplyr' was built under R version 3.6.3## Warning: package 'forcats' was built under R version 3.6.3## -- Conflicts ------------------- tidyverse_conflicts() --
## x dplyr::collapse() masks nlme::collapse()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sandwich)## Warning: package 'sandwich' was built under R version 3.6.3library(lmtest)## Warning: package 'lmtest' was built under R version 3.6.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(modelr)## Warning: package 'modelr' was built under R version 3.6.3library(broom)## Warning: package 'broom' was built under R version 3.6.3##
## Attaching package: 'broom'## The following object is masked from 'package:modelr':
##
## bootstraplmtest
bptest(bb_lm2)##
## studentized Breusch-Pagan test
##
## data: bb_lm2
## BP = 3.8947, df = 4, p-value = 0.4204errores estandar con el paquete sandwich
bb_lm2 %>%
vcovHC() %>%
diag() %>%
sqrt()## (Intercept) energy loudness speechiness valence
## 326.34683 242.31984 18.54877 243.32310 108.57048coeftest(bb_lm2, vcov = vcovHC(bb_lm2, type = "HC0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 726.438 290.950 2.4968 0.016985 *
## energy -533.615 215.641 -2.4745 0.017919 *
## loudness 39.631 16.397 2.4170 0.020559 *
## speechiness 553.657 203.070 2.7264 0.009632 **
## valence -153.147 96.195 -1.5920 0.119661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo lineal con errores robustos
conejo_robusto<-lmrob(Rep ~ energy+loudness+speechiness+valence,correlation =corLin(form=~tempo),data=conejomalo)