#Datos
data("Journals",package = "AER")
#Valida datos
dim(Journals)
## [1] 180 10
#Se ajusta el nombre
journals<-Journals[,c("subs","price")]
#validacion que se creara correctamente
dim(journals)
## [1] 180 2
#comienza el codigo del libro
journals$citerprice<-Journals$price/Journals$citations
summary(journals)
## subs price citerprice
## Min. : 2.0 Min. : 20.0 Min. : 0.005223
## 1st Qu.: 52.0 1st Qu.: 134.5 1st Qu.: 0.464495
## Median : 122.5 Median : 282.0 Median : 1.320513
## Mean : 196.9 Mean : 417.7 Mean : 2.548455
## 3rd Qu.: 268.2 3rd Qu.: 540.8 3rd Qu.: 3.440171
## Max. :1098.0 Max. :2120.0 Max. :24.459460
#continuamos con el codigo del libro pag. 57
plot(log(subs)~log(citerprice),data = journals)
jour_lm <- lm(log(subs)~log(citerprice),data = journals)
abline(jour_lm)
#continua codigo / el texto recomienda usar lm y conocer sus componentes
class(jour_lm)
## [1] "lm"
names(jour_lm)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
#obtenemos los datos de jour_lm
summary(jour_lm)
##
## Call:
## lm(formula = log(subs) ~ log(citerprice), data = journals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.72478 -0.53609 0.03721 0.46619 1.84808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.76621 0.05591 85.25 <2e-16 ***
## log(citerprice) -0.53305 0.03561 -14.97 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7497 on 178 degrees of freedom
## Multiple R-squared: 0.5573, Adjusted R-squared: 0.5548
## F-statistic: 224 on 1 and 178 DF, p-value: < 2.2e-16
#ajustamos lm y obtenemos de nuevo summary
jour_slm <- summary(jour_lm)
class(jour_slm)
## [1] "summary.lm"
names(jour_slm)
## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"
#obtenemos coeficientes del nuevo lm
jour_slm$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7662121 0.05590908 85.24934 2.953913e-146
## log(citerprice) -0.5330535 0.03561320 -14.96786 2.563943e-33
#Analsis de varianza
anova(jour_lm)
## Analysis of Variance Table
##
## Response: log(subs)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(citerprice) 1 125.93 125.934 224.04 < 2.2e-16 ***
## Residuals 178 100.06 0.562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Estimaciones por intervalos
coef(jour_lm)
## (Intercept) log(citerprice)
## 4.7662121 -0.5330535
#se agrega el intervalo
confint(jour_lm,level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 4.6558822 4.8765420
## log(citerprice) -0.6033319 -0.4627751
#Prediccion
predict(jour_lm,newdata = data.frame(citerprice=2.11),interval = "confidence")
## fit lwr upr
## 1 4.368188 4.247485 4.48889
predict(jour_lm,newdata = data.frame(citerprice=2.11),interval = "prediction")
## fit lwr upr
## 1 4.368188 2.883746 5.852629
#Uso de prediccion para visualizar los datos
lciteprice <- seq(from=-6,to=4,by=0.25)
jour_pred <- predict(jour_lm, interval = "prediction",newdata = data.frame(citerprice=exp(lciteprice)))
plot(log(subs)~log(citerprice),data = journals)
lines(jour_pred[,1]~lciteprice,col=1)
lines(jour_pred[,2]~lciteprice,col=1,lty=2)
lines(jour_pred[,3]~lciteprice,col=1,lty=2)
#Graficos para lm
par(mfrow=c(2,2))
plot(jour_lm)
par(mfrow=c(1,1))
#aqui se instala el paquete car que y e agrega pues no lo tenia y sin este la funcion no existe
#Investigando ya que la funcion no daba al parecer cambio a linearHypothesis
linearHypothesis(jour_lm,"log(citerprice)=-0.5")
##
## Linear hypothesis test:
## log(citerprice) = - 0.5
##
## Model 1: restricted model
## Model 2: log(subs) ~ log(citerprice)
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 179 100.54
## 2 178 100.06 1 0.48421 0.8614 0.3546
#Regresion lineal multiple, necesito agregar la base de datos
data("CPS1988",package = "AER")
summary(CPS1988)
## wage education experience ethnicity smsa
## Min. : 50.05 Min. : 0.00 Min. :-4.0 cauc:25923 no : 7223
## 1st Qu.: 308.64 1st Qu.:12.00 1st Qu.: 8.0 afam: 2232 yes:20932
## Median : 522.32 Median :12.00 Median :16.0
## Mean : 603.73 Mean :13.07 Mean :18.2
## 3rd Qu.: 783.48 3rd Qu.:15.00 3rd Qu.:27.0
## Max. :18777.20 Max. :18.00 Max. :63.0
## region parttime
## northeast:6441 no :25631
## midwest :6863 yes: 2524
## south :8760
## west :6091
##
##
#se crea la variable para el ajustado
cps_lm<- lm(log(wage)~experience+I(experience^2)+education+ethnicity,data = CPS1988)
#revisamos los datos de la nueva variable
summary(cps_lm)
##
## Call:
## lm(formula = log(wage) ~ experience + I(experience^2) + education +
## ethnicity, data = CPS1988)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9428 -0.3162 0.0580 0.3756 4.3830
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.321e+00 1.917e-02 225.38 <2e-16 ***
## experience 7.747e-02 8.800e-04 88.03 <2e-16 ***
## I(experience^2) -1.316e-03 1.899e-05 -69.31 <2e-16 ***
## education 8.567e-02 1.272e-03 67.34 <2e-16 ***
## ethnicityafam -2.434e-01 1.292e-02 -18.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5839 on 28150 degrees of freedom
## Multiple R-squared: 0.3347, Adjusted R-squared: 0.3346
## F-statistic: 3541 on 4 and 28150 DF, p-value: < 2.2e-16
#Comparacion de modelos
cps_noeth<- lm(log(wage)~experience+I(experience^2)+education,data = CPS1988)
anova(cps_noeth,cps_lm)
## Analysis of Variance Table
##
## Model 1: log(wage) ~ experience + I(experience^2) + education
## Model 2: log(wage) ~ experience + I(experience^2) + education + ethnicity
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28151 9719.6
## 2 28150 9598.6 1 121.02 354.91 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#esto demuestra como es que la etnia llega a ser significativo en nivel razonable
anova(cps_lm)
## Analysis of Variance Table
##
## Response: log(wage)
## Df Sum Sq Mean Sq F value Pr(>F)
## experience 1 839.5 839.52 2462.06 < 2.2e-16 ***
## I(experience^2) 1 2249.5 2249.49 6597.10 < 2.2e-16 ***
## education 1 1619.7 1619.69 4750.07 < 2.2e-16 ***
## ethnicity 1 121.0 121.02 354.91 < 2.2e-16 ***
## Residuals 28150 9598.6 0.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cps_noeth<- update(cps_lm,formula. = .~.-ethnicity)
waldtest(cps_lm,.~.-ethnicity)
## Wald test
##
## Model 1: log(wage) ~ experience + I(experience^2) + education + ethnicity
## Model 2: log(wage) ~ experience + I(experience^2) + education
## Res.Df Df F Pr(>F)
## 1 28150
## 2 28151 -1 354.91 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#modelos parcialmente lineales
#falta agregar una libreria para que apliquen algunas funciones
cps_plm<-lm(log(wage)~bs(experience,df=5)+education+ethnicity,data = CPS1988)
#ajuste para que este entre 3 y 10
cps_bs<-lapply(3:10, function(i)lm(log(wage)~bs(experience,df=i)+education+ethnicity,data = CPS1988))
structure(sapply(cps_bs,AIC,k=log(nrow(CPS1988))),.Names=3:10)
## 3 4 5 6 7 8 9 10
## 49205.17 48835.97 48794.23 48794.83 48801.41 48797.30 48799.28 48802.49
#comparativa de modelo parcialmente lineal y el clasico
cps<-data.frame(experience=-2.6,education=with(CPS1988,mean(education[ethnicity=="cauc"])),ethnicity="cauc")
cps$yhat1<-predict(cps_lm,newdata = cps)
cps$yhat2<-predict(cps_lm,newdata = cps)
plot(log(wage)~jitter(experience,factor = 3),pch=19,col=rgb(0.5,0.5,0.5,alpha = 0.02),data = CPS1988)
lines(yhat1~experience,data = cps,lty=2)
lines(yhat2~experience,data = cps)
legend("topleft",c("quadratic","spline"),lty=c(2,1),bty="n")
#Factores, Interacciones y Ponderaciones
#Interacciones
cps_int<-lm(log(wage)~experience+I(experience^2)+education*ethnicity,data = CPS1988)
coeftest(cps_int)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3131e+00 1.9590e-02 220.1703 < 2e-16 ***
## experience 7.7520e-02 8.8028e-04 88.0625 < 2e-16 ***
## I(experience^2) -1.3179e-03 1.9006e-05 -69.3388 < 2e-16 ***
## education 8.6312e-02 1.3089e-03 65.9437 < 2e-16 ***
## ethnicityafam -1.2389e-01 5.9026e-02 -2.0989 0.03584 *
## education:ethnicityafam -9.6481e-03 4.6510e-03 -2.0744 0.03805 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#regresiones separadas para cada nivel
cps_sep<-lm(log(wage)~ethnicity/(experience+I(experience^2)+education)-1,data = CPS1988)
cps_sep_cf<-matrix(coef(cps_sep),nrow = 2)
rownames(cps_sep_cf)<-levels(CPS1988$ethnicity)
colnames(cps_sep_cf)<-names(coef(cps_lm))[1:4]
cps_sep_cf
## (Intercept) experience I(experience^2) education
## cauc 4.310196 0.07923367 -0.0013596710 0.08575089
## afam 4.159317 0.06189576 -0.0009414879 0.08654232
#Comparacion del nuevo modelo con el primero
anova(cps_sep,cps_lm)
## Analysis of Variance Table
##
## Model 1: log(wage) ~ ethnicity/(experience + I(experience^2) + education) -
## 1
## Model 2: log(wage) ~ experience + I(experience^2) + education + ethnicity
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28147 9581.8
## 2 28150 9598.6 -3 -16.814 16.464 1.099e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Cambio categoria
CPS1988$region <- relevel(CPS1988$region, ref = "south")
cps_region <- lm(log(wage) ~ ethnicity + education + experience + I(experience^2) + region, data = CPS1988)
coef(cps_region)
## (Intercept) ethnicityafam education experience I(experience^2)
## 4.283606335 -0.225678877 0.084672493 0.077656452 -0.001322942
## regionnortheast regionmidwest regionwest
## 0.131920488 0.043789477 0.040326813
#minimos cuadrados ponderados
jour_wls1 <- lm(log(subs) ~ log(citerprice), data = journals,weights = 1/citerprice^2)
jour_wls2 <- lm(log(subs) ~ log(citerprice), data = journals, weights = 1/citerprice)
plot(log(subs) ~ log(citerprice), data = journals)
abline(jour_lm)
abline(jour_wls1, lwd = 2, lty = 2)
abline(jour_wls2, lwd = 2, lty = 3)
legend("bottomleft", c("OLS", "WLS1", "WLS2"),lty = 1:3, lwd = 2, bty = "n")
auxreg <- lm(log(residuals(jour_lm)^2) ~ log(citerprice),data = journals)
jour_fgls1 <- lm(log(subs) ~ log(citerprice),weights = 1/exp(fitted(auxreg)), data = journals)
gamma2i<-coef(auxreg)[2]
gamma2<- 0
while(abs((gamma2i - gamma2)/gamma2)>1e-7){gamma2<-gamma2i
fglsi<-lm(log(subs)~log(citerprice),data = journals,weights = 1/citerprice^gamma2)
gamma2i<-coef(lm(log(residuals(fglsi)^2)~log(citerprice),data = journals))[2]
}
jour_fgls2<-lm(log(subs)~log(citerprice),data = journals,weights = 1/citerprice^gamma2)
coef(jour_fgls2)
## (Intercept) log(citerprice)
## 4.7758234 -0.5007952
#regresion lineal en series de tiempo
data("USMacroG")
plot(USMacroG[,c("dpi","consumption")],lty=c(3,1),plot.type = "single",ylab="")
legend("topleft",legend=c("income","consumption"),lty=c(3,1),bty="n")
#siguiente ejercicio se instla paquete
data("USMacroG",package = "AER")
cons_lm1 <- dynlm(consumption~dpi+L(dpi),data = USMacroG)
cons_lm2 <- dynlm(consumption~dpi+L(consumption),data = USMacroG)
summary(cons_lm1)
##
## Time series regression with "ts" data:
## Start = 1950(2), End = 2000(4)
##
## Call:
## dynlm(formula = consumption ~ dpi + L(dpi), data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.02 -56.68 1.58 49.91 323.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -81.07959 14.50814 -5.589 7.43e-08 ***
## dpi 0.89117 0.20625 4.321 2.45e-05 ***
## L(dpi) 0.03091 0.20754 0.149 0.882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.58 on 200 degrees of freedom
## Multiple R-squared: 0.9964, Adjusted R-squared: 0.9964
## F-statistic: 2.785e+04 on 2 and 200 DF, p-value: < 2.2e-16
plot(merge(as.zoo(USMacroG[,"consumption"]),fitted(cons_lm1),fitted(cons_lm2),0,residuals(cons_lm1),
residuals(cons_lm2)),screens = rep(1:2,c(3,3)),lty = rep(1:3,2),ylab=c("Fitted values","Residuals"),xlab = "Time",main = "")
legend(0.05,0.95,c("observed","cons_lm1","cons_lm2"),lty = 1:3,bty = "n")
#Encompassing test, agregar libreria
cons_lmE <- dynlm(consumption~dpi+L(dpi)+L(consumption),data = USMacroG)
anova(cons_lm1,cons_lmE,cons_lm2)
## Analysis of Variance Table
##
## Model 1: consumption ~ dpi + L(dpi)
## Model 2: consumption ~ dpi + L(dpi) + L(consumption)
## Model 3: consumption ~ dpi + L(consumption)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 200 1534001
## 2 199 73550 1 1460451 3951.448 < 2.2e-16 ***
## 3 200 92644 -1 -19094 51.661 1.299e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
encomptest(cons_lm1,cons_lm2)
## Encompassing test
##
## Model 1: consumption ~ dpi + L(dpi)
## Model 2: consumption ~ dpi + L(consumption)
## Model E: consumption ~ dpi + L(dpi) + L(consumption)
## Res.Df Df F Pr(>F)
## M1 vs. ME 199 -1 3951.448 < 2.2e-16 ***
## M2 vs. ME 199 -1 51.661 1.299e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#regresion lineal con datos de panel
data("Grunfeld",package = "AER")
gr <- subset(Grunfeld,firm%in%c("General Electric","General Motors","IBM"))
# en la consola recomienda cambiar de plm.data a pdata.frame
pgr <- pdata.frame(gr,index=c("firm","year"))
gr_pool <- plm(invest~value+capital,data = pgr,model = "pooling")
gr_fe <- plm(invest~value+capital,data = pgr,model = "within")
summary(gr_fe)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = invest ~ value + capital, data = pgr, model = "within")
##
## Balanced Panel: n = 3, T = 20, N = 60
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -167.3305 -26.1407 2.0878 26.8442 201.6813
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.104914 0.016331 6.4242 3.296e-08 ***
## capital 0.345298 0.024392 14.1564 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1888900
## Residual Sum of Squares: 243980
## R-Squared: 0.87084
## Adj. R-Squared: 0.86144
## F-statistic: 185.407 on 2 and 55 DF, p-value: < 2.22e-16
pFtest(gr_fe,gr_pool)
##
## F test for individual effects
##
## data: invest ~ value + capital
## F = 56.825, df1 = 2, df2 = 55, p-value = 4.148e-14
## alternative hypothesis: significant effects
gr_re <- plm(invest ~ value + capital, data = pgr,model = "random", random.method = "walhus")
summary(gr_re)
## Oneway (individual) effect Random Effect Model
## (Wallace-Hussain's transformation)
##
## Call:
## plm(formula = invest ~ value + capital, data = pgr, model = "random",
## random.method = "walhus")
##
## Balanced Panel: n = 3, T = 20, N = 60
##
## Effects:
## var std.dev share
## idiosyncratic 4389.31 66.25 0.352
## individual 8079.74 89.89 0.648
## theta: 0.8374
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -187.3987 -32.9206 6.9595 31.4322 210.2006
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -109.976572 61.701384 -1.7824 0.07468 .
## value 0.104280 0.014996 6.9539 3.553e-12 ***
## capital 0.344784 0.024520 14.0613 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1988300
## Residual Sum of Squares: 257520
## R-Squared: 0.87048
## Adj. R-Squared: 0.86594
## Chisq: 383.089 on 2 DF, p-value: < 2.22e-16
plmtest(gr_pool)
##
## Lagrange Multiplier Test - (Honda)
##
## data: invest ~ value + capital
## normal = 15.47, p-value < 2.2e-16
## alternative hypothesis: significant effects
phtest(gr_re,gr_fe)
##
## Hausman Test
##
## data: invest ~ value + capital
## chisq = 0.04038, df = 2, p-value = 0.98
## alternative hypothesis: one model is inconsistent
#modelo lineal dinamico
data("EmplUK",package = "plm")
form<- log(emp)~log(wage)+log(capital)+log(output)
empl_lab <- pgmm(dynformula(form,list(2,1,0,1)),data = EmplUK,index = c("firm","year")
,effect = "twoways",model = "twosteps"
,gmm.inst=~log(emp),lag.gmm=list(c(2,99)))
## Warning in dynformula(form, list(2, 1, 0, 1)): use of 'dynformula()' is
## deprecated, use a multi-part formula instead
## Warning in dynformula(gmm.inst, lag.form = lag.gmm): use of 'dynformula()' is
## deprecated, use a multi-part formula instead
summary(empl_lab)
## Twoways effects Two-steps model Difference GMM
##
## Call:
## pgmm(formula = dynformula(form, list(2, 1, 0, 1)), data = EmplUK,
## effect = "twoways", model = "twosteps", index = c("firm",
## "year"), gmm.inst = ~log(emp), lag.gmm = list(c(2, 99)))
##
## Unbalanced Panel: n = 140, T = 7-9, N = 1031
##
## Number of Observations Used: 611
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.6190677 -0.0255683 0.0000000 -0.0001339 0.0332013 0.6410272
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## lag(log(emp), c(1, 2))1 0.474151 0.185398 2.5575 0.0105437 *
## lag(log(emp), c(1, 2))2 -0.052967 0.051749 -1.0235 0.3060506
## log(wage) -0.513205 0.145565 -3.5256 0.0004225 ***
## lag(log(wage), 1) 0.224640 0.141950 1.5825 0.1135279
## log(capital) 0.292723 0.062627 4.6741 2.953e-06 ***
## log(output) 0.609775 0.156263 3.9022 9.530e-05 ***
## lag(log(output), 1) -0.446373 0.217302 -2.0542 0.0399605 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sargan test: chisq(25) = 30.11247 (p-value = 0.22011)
## Autocorrelation test (1): normal = -1.53845 (p-value = 0.12394)
## Autocorrelation test (2): normal = -0.2796829 (p-value = 0.77972)
## Wald test for coefficients: chisq(7) = 142.0353 (p-value = < 2.22e-16)
## Wald test for time dummies: chisq(6) = 16.97046 (p-value = 0.0093924)
#Sistema de ecuaciones lineales
gr2 <- subset(Grunfeld,firm%in%c("Chrysler","IBM"))
pgr2 <- pdata.frame(gr2,c("firm","year"))
gr_sur <- systemfit(invest~value+capital,method = "SUR",data = pgr2)
summary(gr_sur,residCov=FALSE,equations=FALSE)
##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 40 34 4113.64 11022 0.928939 0.927094
##
## N DF SSR MSE RMSE R2 Adj R2
## Chrysler 20 17 3001.64 176.5672 13.28786 0.913457 0.903276
## IBM 20 17 1112.00 65.4117 8.08775 0.952079 0.946441
##
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Chrysler_(Intercept) -5.7031286 13.2773843 -0.42954 0.67292712
## Chrysler_value 0.0779871 0.0195817 3.98266 0.00096271 ***
## Chrysler_capital 0.3114785 0.0286957 10.85455 4.5944e-09 ***
## IBM_(Intercept) -8.0908189 4.5216484 -1.78935 0.09138881 .
## IBM_value 0.1272417 0.0306021 4.15794 0.00065881 ***
## IBM_capital 0.0966341 0.0983302 0.98275 0.33951047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1