Mahalanobis

Primero transformé las variables originales con Box-Cox para aplicar distancias de mahalanobis en la identificación de outliers.

temp <- base[17:47+10]
temp[22:31] <- temp[22:31]+1


for(i in 1:31){
  bc <- boxcox(unlist(temp[i]) ~ 1,data=temp, plotit = F)
  lamda<-bc$x[which.max(bc$y)]
  temp[i]<-ifelse(lamda==0,log(temp[i]),(temp[i]^lamda-1)/lamda)
}

temp$mahalanobis <- mahalanobis(temp, colMeans(temp), cov(temp))
temp$pvalue <- pchisq(temp$mahalanobis, df=ncol(temp)-1, lower.tail=FALSE)

base <- base[which(temp$pvalue>0.05),]

rm(list = c("temp", "bc", "i","lamda"))

A partir de las distancias de mahalanobis se detectaron 189 casos cuyas distancias multivariadas se alejaron significativamente del resto de datos, quedando un total de 941 casos.

Cualitativo

Se realizó un análisis cualitativo de los puntajes de audit para evaluar la concordancia entre respuestas (por ejemplo, si nunca consumió alcohol, no puede presentar CEEA más de una vez al mes).

base <- base[!(base$audit01==0&
       (base$audit02>1|
          base$audit02>1|
          base$audit03>1|
          base$audit04>1|
          base$audit05>1|
          base$audit06>1|
          base$audit08>1)),]

A partir de este análisis se eliminaron 27 casos con respuestas incongruentes, resultando en una base final de 914 casos.

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_1~AT+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad
## Model 2: alcohol_1 ~ AT + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     250.18                     
## 2       911     250.15  1 0.027723   0.8678

# ajuste<-glm(alcohol_1~uni+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~ben+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo, pero AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_1~AP+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad
## Model 2: alcohol_1 ~ AP + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     250.18                     
## 2       911     249.88  1  0.29739   0.5855

# 
# ajuste<-glm(alcohol_1~pod+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~log+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo, pero AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_1~AC+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad
## Model 2: alcohol_1 ~ AC + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       912     250.18                        
## 2       911     240.58  1   9.5949 0.001951 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_1 ~ AC + edad, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -20.0257     3.0759  -6.510 7.49e-11 ***
## AC            0.8975     0.3110   2.886   0.0039 ** 
## edad          0.7431     0.1409   5.274 1.33e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 240.58  on 911  degrees of freedom
## AIC: 246.58
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

## Waiting for profiling to be done...

# ajuste<-glm(alcohol_1~hed+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~est+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~adi+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo y C fue significativo en la predicción de consumo de alcohol (\(OR\) = 2.45; \(IC_{95}\) = [1.37; 4.66]). De acuerdo a este resultado por cada unidad que aumenta la apertura al cambio, el riesgo de presentar una frecuencia de consumo elevada aumenta en un 145.34%.

htest <- hoslem.test(base$alcohol_1, ajuste$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_1, ajuste$fitted.values
## X-squared = 2.7932e-19, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_1~AC+edad, data=base, family=binomial, subset=train)

clase<-base$alcohol_1[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_1[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_1[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 225  43
##   1   2   5

Pred.: Conservación (C)

ajuste<-glm(alcohol_1~C+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad
## Model 2: alcohol_1 ~ C + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       912     250.18                        
## 2       911     240.07  1   10.107 0.001477 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_1 ~ C + edad, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -12.6218     2.4445  -5.163 2.43e-07 ***
## C            -0.7283     0.2376  -3.065  0.00218 ** 
## edad          0.7335     0.1401   5.235 1.65e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 240.07  on 911  degrees of freedom
## AIC: 246.07
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

## Waiting for profiling to be done...

# ajuste<-glm(alcohol_1~con+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~tra+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_1~seg+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# summary(glm(alcohol_1~seg+tra+con+edad+sexo, data=base, family=binomial))

Se cumplen los supuestos del modelo y AC fue significativo en la predicción de consumo de alcohol (\(OR\) = 0.48; \(IC_{95}\) = [0.3; 0.76]). De acuerdo a este resultado por cada unidad que aumenta la conservación, el riesgo de presentar una frecuencia de consumo elevada disminuye en un 51.73%.

htest <- hoslem.test(base$alcohol_1, ajuste$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_1, ajuste$fitted.values
## X-squared = 1.3016e-19, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_1~edad+C, data=base, family=binomial, subset=train)

clase<-base$alcohol_1[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_1[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_1[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 198  70
##   1   3   4

Modelo final (Audit 1)

Se tomó como punto de base el modelo con edad y conservación por ser el que obtuvo una menor deviance y se comparó con el modelo que además tenía apertura al cambio.

mod_C<-glm(alcohol_1~edad+C, data=base, family=binomial)
mod_total <- glm(alcohol_1~edad+C+AC, data=base, family=binomial)
# plot(simulateResiduals(mod_total))
anova(mod_C, mod_total, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad + C
## Model 2: alcohol_1 ~ edad + C + AC
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       911     240.07                          
## 2       910     225.92  1   14.145 0.0001693 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod_total)

## 
## Call:
## glm(formula = alcohol_1 ~ edad + C + AC, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -18.5657     3.2155  -5.774 7.75e-09 ***
## edad          0.7777     0.1484   5.239 1.61e-07 ***
## C            -0.8568     0.2352  -3.644 0.000269 ***
## AC            1.1232     0.3270   3.435 0.000593 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 225.92  on 910  degrees of freedom
## AIC: 233.92
## 
## Number of Fisher Scoring iterations: 7

OR_edad <- exp(mod_total$coefficients[2])
OR_C <- exp(mod_total$coefficients[3])
OR_AC <- exp(mod_total$coefficients[4])
IC_OR_edad<-exp(confint(mod_total))[2,]

## Waiting for profiling to be done...

IC_OR_C<-exp(confint(mod_total))[3,]

## Waiting for profiling to be done...

IC_OR_AC<-exp(confint(mod_total))[4,]

## Waiting for profiling to be done...

R2 <- NagelkerkeR2(mod_total)$R2

Se cumplen los supuestos del modelo y edad, conservación y apertura al cambio fueron significativos en la predicción de consumo de alcohol (\(OR_{edad}\) = 2.18, \(IC_{95}\) = [1.65; 2.97]; \(OR_{Conservacion}\) = 0.42, \(IC_{95}\) = [0.26; 0.66]; \(OR_{AperturaAlCambio}\) = 3.07, \(IC_{95}\) = [1.67; 6.05]). De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar una frecuencia alta de consumo es 117.65% mayor por cada año que aumenta la edad, 207.45% mayor por cada unidad que aumenta la apertura al cambio y 57.55% menor por cada unidad que aumenta la conservación.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.231) se estima que el modelo propuesto explica el 23.05% de la variabilidad de Audit 1.

htest <- hoslem.test(base$alcohol_1, mod_total$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_1, mod_total$fitted.values
## X-squared = 4.3918e-17, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_1~edad+C+AC, data=base, family=binomial, subset=train)

clase<-base$alcohol_1[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_1[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_1[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 193  75
##   1   2   5

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_2~AT+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AT + edad + sexo
##   Resid. Df Resid. Dev Df  Deviance Pr(>Chi)
## 1       911     431.34                      
## 2       910     431.34  1 2.488e-05    0.996

# ajuste<-glm(alcohol_2~uni+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~ben+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo, pero AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_2~AP+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AP + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     431.34                     
## 2       910     431.29  1 0.047318   0.8278

# ajuste<-glm(alcohol_2~pod+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~log+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo, pero AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_2~AC+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AC + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     431.34                     
## 2       910     431.05  1  0.29079   0.5897

# ajuste<-glm(alcohol_2~adi+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~hed+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~est+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")

Se cumplen los supuestos del modelo, pero AC no resultó significativo en la predicción de consumo de alcohol.

Pred.: Conservación (C)

ajuste<-glm(alcohol_2~C+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ C + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       911     431.34                       
## 2       910     426.01  1   5.3225  0.02105 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_2 ~ C + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -13.40995    1.72878  -7.757 8.70e-15 ***
## C            -0.35935    0.15800  -2.274  0.02295 *  
## edad          0.77434    0.09981   7.758 8.64e-15 ***
## sexom        -0.73349    0.27171  -2.700  0.00694 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 513.88  on 913  degrees of freedom
## Residual deviance: 426.01  on 910  degrees of freedom
## AIC: 434.01
## 
## Number of Fisher Scoring iterations: 6

OR_C <- exp(ajuste$coefficients[2])
OR_edad <- exp(ajuste$coefficients[3])
OR_sexo <- exp(ajuste$coefficients[4])
IC_OR_C<-exp(confint(ajuste))[2,]

## Waiting for profiling to be done...

IC_OR_edad<-exp(confint(ajuste))[3,]

## Waiting for profiling to be done...

IC_OR_sexo<-exp(confint(ajuste))[4,]

## Waiting for profiling to be done...

R2 <- NagelkerkeR2(ajuste)$R2


# ajuste<-glm(alcohol_2~con+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~tra+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# ajuste<-glm(alcohol_2~seg+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
# anova(mod_edad, ajuste, test = "Chisq")
# 
# summary(glm(alcohol_2~seg+tra+con+edad+sexo, data=base, family=binomial))

Se cumplen los supuestos del modelo y la edad, el sexo y la conservación fueron significativos en la predicción de consumo de alcohol (\(OR_{edad}\) = 2.17, \(IC_{95}\) = [1.8; 2.66]; \(OR_{sexo}\) = 0.48, \(IC_{95}\) = [0.28; 0.81]; \(OR_{Conservación}\) = 0.7, \(IC_{95}\) = [0.51; 0.95]). De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar CEEA es 116.92% mayor por cada año que aumenta la edad, 30.19% menor por cada unidad que aumenta la conservación y 51.98% en varones que en mujeres.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.213) se estima que el modelo propuesto explica el 21.31% de la variabilidad de Audit 2.

htest <- hoslem.test(base$alcohol_2, ajuste$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_2, ajuste$fitted.values
## X-squared = 2.5847e-17, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_2~C+edad+sexo, data=base, family=binomial, subset=train)

clase<-base$alcohol_2[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_2[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_2[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 181  71
##   1   6  17

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_3~AT+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ AT + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     240.35                     
## 2       911     238.87  1   1.4873   0.2226

ajuste<-glm(alcohol_3~uni+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ uni + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     240.35                     
## 2       911     239.42  1  0.92964    0.335

ajuste<-glm(alcohol_3~ben+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ ben + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     240.35                     
## 2       911     238.87  1   1.4846   0.2231

ajuste<-glm(alcohol_3~seg+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ seg + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       912     240.35                     
## 2       911     238.32  1   2.0337   0.1538

summary(glm(alcohol_3~seg+tra+con+edad, data=base, family=binomial))

## 
## Call:
## glm(formula = alcohol_3 ~ seg + tra + con + edad, family = binomial, 
##     data = base)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -16.67703    2.74207  -6.082 1.19e-09 ***
## seg          -0.15668    0.16586  -0.945    0.345    
## tra          -0.20266    0.19023  -1.065    0.287    
## con          -0.08696    0.16803  -0.518    0.605    
## edad          0.91452    0.15407   5.936 2.93e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 236.47  on 909  degrees of freedom
## AIC: 246.47
## 
## Number of Fisher Scoring iterations: 7

Hay un levísimo incumplimiento de supuestos, pero igual AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_3~AP+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ AP + edad
##   Resid. Df Resid. Dev Df  Deviance Pr(>Chi)
## 1       912     240.35                      
## 2       911     240.35  1 0.0052026   0.9425

Hay un levísimo incumplimiento de supuestos, pero igual AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_3~AC+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ AC + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       912     240.35                        
## 2       911     233.60  1   6.7481 0.009385 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_3 ~ AC + edad, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.5057     3.2920  -6.836 8.12e-12 ***
## AC            0.7375     0.2996   2.461   0.0139 *  
## edad          0.9374     0.1559   6.014 1.81e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 233.60  on 911  degrees of freedom
## AIC: 239.6
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

## Waiting for profiling to be done...

Se cumplen los supuestos del modelo y AC fue significativo en la predicción de consumo de alcohol (\(OR\) = 2.09; \(IC_{95}\) = [1.19; 3.87]). De acuerdo a este resultado por cada unidad que aumenta la apertura al cambio, el riesgo de presentar una frecuencia de CEEA elevada aumenta en un 109.06%.

htest <- hoslem.test(base$alcohol_3, ajuste$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_3, ajuste$fitted.values
## X-squared = 1.4792e-16, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_3~AC+edad, data=base, family=binomial, subset=train)

clase<-base$alcohol_3[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_3[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_3[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 203  63
##   1   3   6

Pred.: Conservación (C)

ajuste<-glm(alcohol_3~C+edad, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ C + edad
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       912     240.35                       
## 2       911     236.64  1   3.7085  0.05413 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_3 ~ C + edad, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -16.7865     2.7219  -6.167 6.95e-10 ***
## C            -0.4328     0.2290  -1.889   0.0588 .  
## edad          0.9156     0.1535   5.967 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 236.64  on 911  degrees of freedom
## AIC: 242.64
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

## Waiting for profiling to be done...

Se cumplen los supuestos del modelo y C tuvo una significación marginal en la predicción de consumo de alcohol (\(OR\) = 0.65; \(IC_{95}\) = [0.41; 1.01]). De acuerdo a este resultado por cada unidad que aumenta la conservación, el riesgo de presentar una frecuencia de consumo elevada disminuye en un 35.13%.

htest <- hoslem.test(base$alcohol_3, ajuste$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_3, ajuste$fitted.values
## X-squared = 3.2757e-17, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_3~C+edad, data=base, family=binomial, subset=train)

clase<-base$alcohol_3[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_3[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_3[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 219  47
##   1   4   5

Modelo final (Audit 3)

Se tomó como punto de base el modelo con edad y apertura al cambio por ser el que obtuvo una menor deviance y se comparó con el modelo que además tenía conservación.

mod_AC<-glm(alcohol_3~edad+AC, data=base, family=binomial)
mod_total <- glm(alcohol_3~edad+AC+C, data=base, family=binomial)
# plot(simulateResiduals(mod_total))
anova(mod_AC, mod_total, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad + AC
## Model 2: alcohol_3 ~ edad + AC + C
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       911     233.60                       
## 2       910     227.44  1   6.1663  0.01302 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod_total)

## 
## Call:
## glm(formula = alcohol_3 ~ edad + AC + C, family = binomial, data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -21.6098     3.3828  -6.388 1.68e-10 ***
## edad          0.9636     0.1610   5.984 2.18e-09 ***
## AC            0.8842     0.3109   2.844  0.00446 ** 
## C            -0.5480     0.2257  -2.429  0.01516 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 227.44  on 910  degrees of freedom
## AIC: 235.44
## 
## Number of Fisher Scoring iterations: 7

OR_edad <- exp(mod_total$coefficients[2])
OR_AC <- exp(mod_total$coefficients[3])
OR_C <- exp(mod_total$coefficients[4])
IC_OR_edad<-exp(confint(mod_total))[2,]

## Waiting for profiling to be done...

IC_OR_AC<-exp(confint(mod_total))[3,]

## Waiting for profiling to be done...

IC_OR_C<-exp(confint(mod_total))[4,]

## Waiting for profiling to be done...

R2 <- NagelkerkeR2(mod_total)$R2

Se cumplen los supuestos del modelo y edad, conservación y apertura al cambio fueron significativos en la predicción de consumo de alcohol (\(OR_{edad}\) = 2.62, \(IC_{95}\) = [1.95; 3.68]; \(OR_{Conservacion}\) = 0.58, \(IC_{95}\) = [0.37; 0.89]; \(OR_{AperturaAlCambio}\) = 2.42, \(IC_{95}\) = [1.35; 4.59]). De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar una frecuencia alta de CEEA es 162.12% mayor por cada año que aumenta la edad, 142.09% mayor por cada unidad que aumenta la apertura al cambio y 42.19% menor por cada unidad que aumenta la conservación.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.245) se estima que el modelo propuesto explica el 24.5% de la variabilidad de Audit 3.

htest <- hoslem.test(base$alcohol_3, mod_total$fitted.values, g=1)

## Warning in pchisq(chisq, g - 2): NaNs produced

htest

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base$alcohol_3, mod_total$fitted.values
## X-squared = 1.5205e-15, df = -1, p-value = NA

sum(htest$expected<5) # Si hay menores que 5 tengo que rearmar los grupos

## [1] 0

# con el argumento "g" para que sea válido el test (arranqué en 10 y fui bajando)

set.seed(195)
train<-sample(1:nrow(base), nrow(base)*0.7)

ajuste2<-glm(alcohol_3~edad+C+AC, data=base, family=binomial, subset=train)

clase<-base$alcohol_3[-train]
ajustados<-predict(ajuste2,newdata = base[-train,],type="response")
ROC <- rocit(score=ajustados,class=clase)
distancia <- sqrt(ROC$FPR^2 +(ROC$TPR-1)^2)
plot(ROC)

plot(ROC$Cutoff,distancia,pch=20)

minimo <- ROC$Cutoff[which.min(distancia)]
tita<-minimo
ajustados1<-predict(ajuste2,newdata = base[-train,],type="response")
predichos1<-ifelse(ajustados1>tita,1,0)
confusion<-table(base$alcohol_3[-train], predichos1)
tasa_error_te1<-mean(base$alcohol_3[-train]!=predichos1) 

sensibilidad<-confusion[2,2]/sum(confusion[2,]) # capacidad del modelo para detectar los positivos
especificidad<-confusion[1,1]/sum(confusion[1,]) #capacidad del modelo para detectar los negativos

confusion

##    predichos1
##       0   1
##   0 189  77
##   1   2   7

Edad

Como era de esperarse, en todos los modelos para predecir frecuencia de consumo, presencia de CEEA y frecuencia de CEEA la edad fue significativa, con ORs en los modelos finales de 2.18 para la frecuencia de consumo, 2.17 para la presencia de CEEA y 2.62 para la frecuencia de CEEA. De esta forma, la edad fue el único predictor significativo en todos los modelos, siendo que por cada año que aumenta la edad, el riesgo de presentar una frecuencia de consumo elevada, de presencia de CEEA y de frecuencia de CEEA elevada aumentan más del doble.

Sexo

El sexo fue significativo sólo en los cuatro modelos de predicción de la presencia de CEEA, con un OR en el modelo final de 0.48 para los varones. Es decir, las chances de presentar un patrón de consumo de riesgo (CEEA) en varones es un 51.98% menor que en mujeres. En la última EMSE algo de esto aparece, aunque no recuerdo si específicamente con la presencia de CEEA.

Valores

Los valores de autotrascendencia y de autopromoción no fueron predictores significativos en ninguno de los modelos.

Los valores de apertura al cambio permitieron predecir la probabilidad de presentar frecuencia alta de consumo (OR = 3.07) y frecuencia alta de CEEA (OR = 2.42), pero no presencia de CEEA. Que los motivos por los que se comportan como se comportan estén vinculados a la búsqueda de autonomía e independencia, de estimulación, y de disfrute y placer aumenta el riesgo de presentar una frecuencia de consumo y una frecuencia de CEEA elevadas.

Los valores de conservación (i.e., seguridad, tradición y conformidad) predijeron significativamente la probabilidad de frecuencia de consumo (OR = 0.42), la presencia de CEEA (OR = 0.7) y la frecuencia de CEEA (OR = 0.58). Es decir, al tener niveles más altos de conservación, el riesgo de tener una frecuencia de consumo elevada, de presentar CEEA y de tener una frecuencia de CEEA elevada disminuye.

Modelos finales

Conclusión MDS

A pesar de los pequeños inconvenientes, la estructura general respeta la distribución teórica. Esto y el indicador Stress-1 aceptable permiten pasar al siguiente paso (CFA con estrategia de lupa). Dada la configuración obtenida, los valores de hedonismo serán considerados como parte de Apertura al cambio.

Autotrascendencia

mod <- ' universalism =~ un03 + un08 + un19
         benevolence =~ be12 + be18
        '
fit_ml <- cfa(mod, estimator="ML", data=base)
fit_uls <- cfa(mod, estimator="ULS", data=base)


tabla_ajustes_AT <- as.data.frame(rbind(fitmeasures(fit_ml)[c("cfi","rmsea","srmr","rmsea.pvalue")],
                    fitmeasures(fit_uls)[c("cfi","rmsea","srmr","rmsea.pvalue")]))
tabla_ajustes_AT <- cbind(c("ML","ULS"),
                       c(fitmeasures(fit_ml)[c("chisq")]/fitmeasures(fit_ml)[c("df")],
                         fitmeasures(fit_uls)[c("chisq")]/fitmeasures(fit_uls)[c("df")]),
                       tabla_ajustes_AT)
colnames(tabla_ajustes_AT) <- c("Estimador","X2/gl","CFI","RMSEA","SRMR", "PCLOSE")
tabla_ajustes_AT[,2:6] <- round(tabla_ajustes_AT[,2:6],3)

kable(tabla_ajustes_AT,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

Si bien ambos ajustes obtuvieron muy buenos indicadores de bondad de ajuste absolutos (\(SRMR \le 0.08\)), incrementales (\(CFI \ge 0.95\)) y parsimoniosos (\(RMSEA \le 0.05\) y \(PCLOSE \ge 0.05\)) de acuerdo con los criterios de Gana y Broc (2019), el ajuste con ULS superó en bondad de ajuste al realizado con ML en la mayoría de los indicadores.

summary(fit_uls, fit.measures = F, standardized=T)

## lavaan 0.6.15 ended normally after 26 iterations
## 
##   Estimator                                        ULS
##   Optimization method                           NLMINB
##   Number of model parameters                        11
## 
##   Number of observations                           914
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.095
##   Degrees of freedom                                 4
##   P-value (Unknown)                                 NA
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   universalism =~                                                       
##     un03              1.000                               0.667    0.619
##     un08              1.022    0.074   13.805    0.000    0.682    0.582
##     un19              1.025    0.074   13.805    0.000    0.684    0.581
##   benevolence =~                                                        
##     be12              1.000                               0.637    0.728
##     be18              0.739    0.065   11.370    0.000    0.471    0.594
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   universalism ~~                                                       
##     benevolence       0.360    0.024   14.851    0.000    0.847    0.847
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .un03              0.716    0.053   13.425    0.000    0.716    0.617
##    .un08              0.910    0.055   16.693    0.000    0.910    0.662
##    .un19              0.918    0.055   16.783    0.000    0.918    0.662
##    .be12              0.359    0.066    5.432    0.000    0.359    0.470
##    .be18              0.407    0.046    8.942    0.000    0.407    0.648
##     universalism      0.445    0.042   10.638    0.000    1.000    1.000
##     benevolence       0.406    0.057    7.088    0.000    1.000    1.000

plot_AT <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,7)), 
          edge.label.cex = c(rep(1,11)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, 
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

Autopromoción

mod <- ' power =~ po02 + po17
         achievement =~ ac13 + ac04
        '
fit_ml <- cfa(mod, estimator="ML", data=base)
fit_uls <- cfa(mod, estimator="ULS", data=base)


tabla_ajustes_AP <- as.data.frame(rbind(fitmeasures(fit_ml)[c("cfi","rmsea","srmr","rmsea.pvalue")],
                    fitmeasures(fit_uls)[c("cfi","rmsea","srmr","rmsea.pvalue")]))
tabla_ajustes_AP <- cbind(c("ML","ULS"),
                       c(fitmeasures(fit_ml)[c("chisq")]/fitmeasures(fit_ml)[c("df")],
                         fitmeasures(fit_uls)[c("chisq")]/fitmeasures(fit_uls)[c("df")]),
                       tabla_ajustes_AP)
colnames(tabla_ajustes_AP) <- c("Estimador","X2/gl","CFI","RMSEA","SRMR", "PCLOSE")
tabla_ajustes_AP[,2:6] <- round(tabla_ajustes_AP[,2:6],3)

kable(tabla_ajustes_AP,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

Si bien ambos ajustes obtuvieron muy buenos indicadores de bondad de ajuste absolutos (\(SRMR \le 0.08\)), incrementales (\(CFI \ge 0.95\)) y parsimoniosos (\(RMSEA \le 0.05\) y \(PCLOSE \ge 0.05\)) de acuerdo con los criterios de Gana y Broc (2019), el ajuste con ML superó en bondad de ajuste al realizado con ULS en la mayoría de los indicadores.

summary(fit_uls, fit.measures = F, standardized=T)

## lavaan 0.6.15 ended normally after 24 iterations
## 
##   Estimator                                        ULS
##   Optimization method                           NLMINB
##   Number of model parameters                         9
## 
##   Number of observations                           914
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 1.375
##   Degrees of freedom                                 1
##   P-value (Unknown)                                 NA
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   power =~                                                              
##     po02              1.000                               0.891    0.569
##     po17              1.080    0.041   26.233    0.000    0.962    0.617
##   achievement =~                                                        
##     ac13              1.000                               1.095    0.795
##     ac04              0.953    0.036   26.280    0.000    1.044    0.727
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   power ~~                                                              
##     achievement       0.856    0.029   30.009    0.000    0.878    0.878
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .po02              1.662    0.054   30.593    0.000    1.662    0.677
##    .po17              1.502    0.060   24.981    0.000    1.502    0.619
##    .ac13              0.696    0.066   10.513    0.000    0.696    0.367
##    .ac04              0.974    0.062   15.783    0.000    0.974    0.472
##     power             0.794    0.043   18.429    0.000    1.000    1.000
##     achievement       1.199    0.057   20.916    0.000    1.000    1.000

plot_AP <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,4),.75,1), 
          edge.label.cex = c(rep(1,9)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, 
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

Apertura al cambio

mod <- ' selfdirection =~ sd01 + sd11 
         stimulation =~ st06 + st15
         hedonism =~ he10 + he21
        '
fit_ml <- cfa(mod, estimator="ML", data=base)
fit_uls <- cfa(mod, estimator="ULS", data=base)


tabla_ajustes_AC <- as.data.frame(rbind(fitmeasures(fit_ml)[c("cfi","rmsea","srmr","rmsea.pvalue")],
                    fitmeasures(fit_uls)[c("cfi","rmsea","srmr","rmsea.pvalue")]))
tabla_ajustes_AC <- cbind(c("ML","ULS"),
                       c(fitmeasures(fit_ml)[c("chisq")]/fitmeasures(fit_ml)[c("df")],
                         fitmeasures(fit_uls)[c("chisq")]/fitmeasures(fit_uls)[c("df")]),
                       tabla_ajustes_AC)
colnames(tabla_ajustes_AC) <- c("Estimador","X2/gl","CFI","RMSEA","SRMR", "PCLOSE")
tabla_ajustes_AC[,2:6] <- round(tabla_ajustes_AC[,2:6],3)

kable(tabla_ajustes_AC,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

Ambos ajustes obtuvieron buenos indicadores de bondad de ajuste absolutos (SRMR) e incrementales (CFI). En cuanto a los indicadortes parsimoniosos, el ajuste con ULS (al límite de ser buen ajuste) obtuvo mejores niveles de RMSEA que ML (ajuste pobre).

summary(fit_uls, fit.measures = F, standardized=T)

## lavaan 0.6.15 ended normally after 32 iterations
## 
##   Estimator                                        ULS
##   Optimization method                           NLMINB
##   Number of model parameters                        15
## 
##   Number of observations                           914
## 
## Model Test User Model:
##                                                       
##   Test statistic                                41.547
##   Degrees of freedom                                 6
##   P-value (Unknown)                                 NA
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   selfdirection =~                                                      
##     sd01              1.000                               0.670    0.570
##     sd11              0.863    0.049   17.506    0.000    0.578    0.596
##   stimulation =~                                                        
##     st06              1.000                               0.825    0.682
##     st15              1.189    0.064   18.563    0.000    0.981    0.766
##   hedonism =~                                                           
##     he10              1.000                               0.516    0.589
##     he21              1.445    0.110   13.187    0.000    0.746    0.763
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   selfdirection ~~                                                      
##     stimulation       0.494    0.025   19.408    0.000    0.894    0.894
##     hedonism          0.252    0.020   12.485    0.000    0.729    0.729
##   stimulation ~~                                                        
##     hedonism          0.266    0.020   13.240    0.000    0.624    0.624
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .sd01              0.934    0.057   16.449    0.000    0.934    0.675
##    .sd11              0.608    0.048   12.739    0.000    0.608    0.645
##    .st06              0.780    0.057   13.764    0.000    0.780    0.534
##    .st15              0.678    0.073    9.290    0.000    0.678    0.414
##    .he10              0.501    0.045   11.132    0.000    0.501    0.653
##    .he21              0.401    0.072    5.575    0.000    0.401    0.419
##     selfdirection     0.449    0.046    9.732    0.000    1.000    1.000
##     stimulation       0.680    0.046   14.781    0.000    1.000    1.000
##     hedonism          0.267    0.031    8.726    0.000    1.000    1.000

plot_AC <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,9)), 
          edge.label.cex = c(rep(1,15)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, 
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

Conservación

mod <- ' security =~ se05 + se14
         conformity =~ co07 + co16
         tradition =~ tr09 + tr20
        '
fit_ml <- cfa(mod, estimator="ML", data=base)
fit_uls <- cfa(mod, estimator="ULS", data=base)


tabla_ajustes_C <- as.data.frame(rbind(fitmeasures(fit_ml)[c("cfi","rmsea","srmr","rmsea.pvalue")],
                    fitmeasures(fit_uls)[c("cfi","rmsea","srmr","rmsea.pvalue")]))
tabla_ajustes_C <- cbind(c("ML","ULS"),
                       c(fitmeasures(fit_ml)[c("chisq")]/fitmeasures(fit_ml)[c("df")],
                         fitmeasures(fit_uls)[c("chisq")]/fitmeasures(fit_uls)[c("df")]),
                       tabla_ajustes_C)
colnames(tabla_ajustes_C) <- c("Estimador","X2/gl","CFI","RMSEA","SRMR", "PCLOSE")
tabla_ajustes_C[,2:6] <- round(tabla_ajustes_C[,2:6],3)

kable(tabla_ajustes_C,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Indicadores de ajuste según estimador") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

Ambos ajustes obtuvieron buenos indicadores de bondad de ajuste absolutos (SRMR) e incrementales (CFI). En cuanto a los indicadortes parsimoniosos, el ajuste con ML obtuvo mejores niveles de RMSEA que ULS (ajuste pobre).

summary(fit_uls, fit.measures = F, standardized=T)

## lavaan 0.6.15 ended normally after 33 iterations
## 
##   Estimator                                        ULS
##   Optimization method                           NLMINB
##   Number of model parameters                        15
## 
##   Number of observations                           914
## 
## Model Test User Model:
##                                                       
##   Test statistic                                77.613
##   Degrees of freedom                                 6
##   P-value (Unknown)                                 NA
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   security =~                                                           
##     se05              1.000                               0.667    0.503
##     se14              1.294    0.070   18.500    0.000    0.863    0.556
##   conformity =~                                                         
##     co07              1.000                               0.809    0.555
##     co16              1.187    0.055   21.554    0.000    0.961    0.640
##   tradition =~                                                          
##     tr09              1.000                               0.390    0.307
##     tr20              1.644    0.095   17.313    0.000    0.641    0.405
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   security ~~                                                           
##     conformity        0.397    0.022   18.379    0.000    0.736    0.736
##     tradition         0.243    0.017   14.443    0.000    0.935    0.935
##   conformity ~~                                                         
##     tradition         0.348    0.021   16.835    0.000    1.103    1.103
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .se05              1.316    0.048   27.282    0.000    1.316    0.747
##    .se14              1.668    0.067   24.722    0.000    1.668    0.691
##    .co07              1.471    0.053   27.824    0.000    1.471    0.692
##    .co16              1.329    0.067   19.874    0.000    1.329    0.590
##    .tr09              1.456    0.040   36.665    0.000    1.456    0.906
##    .tr20              2.087    0.068   30.709    0.000    2.087    0.836
##     security          0.445    0.035   12.670    0.000    1.000    1.000
##     conformity        0.654    0.041   15.878    0.000    1.000    1.000
##     tradition         0.152    0.022    6.918    0.000    1.000    1.000

plot_C <- semPaths(# Argumentos globales
          fit_uls, what="diagram", whatLabels="std",layout="tree3", rotation = 4, width=50, height=35,exoVar = F,
         # Etiquetas
          label.cex=c(rep(1,9)), 
          edge.label.cex = c(rep(1,15)),
         # Forma de los nodos (latentes)
          nCharNodes = 0, sizeLat = 15,sizeLat2 = 10,border.width=2,
         # Forma de las flechas
          edge.width=2,asize=2,curvePivot=T, 
         # Color
         border.color = "#BFB29E", edge.color = "#BFB29E", edge.label.bg = "#141415", title.color = "#BFB29E", label.color="#BFB29E",bg="transparent", trans=T, vTrans=0
         
         )

Conclusión CFAs

Los cuatro modelos presentaron muy buenos indicadores de bondad de ajuste. Sin embargo, como en Beramendi y Zubieta (2017), algunos valores presentaron una fuerte colinealidad. En este estudio, las asociaciones poder-logro (SE), autodirección-estimulación (OC), seguridad-tradición (C) y conformidad-tradición (C) presentaron cargas factoriales estandarizadas mayores a 0.85, indicador de multicolinealidad. En el caso específico de la asociación entre conformidad y tradición, este coeficiente superó el rango válido [-1; 1], aunque de acuerdo a Deegan (1978) esto….

knitr::include_graphics("zubietaCFA.png")

tabla_ajustes <- as.data.frame(rbind(tabla_ajustes_C[2,],
                                     tabla_ajustes_AT[2,],
                                     tabla_ajustes_AC[2,],
                                     tabla_ajustes_AP[2,]
                                     ))
tabla_ajustes[1] <- c("Conservation",
                      "Self-Transcendence",
                      "Openness to Change",
                      "Self-Enhancement")
colnames(tabla_ajustes)[1] <- "Modelo"

kable(tabla_ajustes,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Indicadores de ajuste para cada dimensión") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

knitr::include_graphics("zubietaCFA_p.png")

par(mfrow = c(2, 2))
plot(plot_AC)
plot(plot_C)
plot(plot_AT)
plot(plot_AP)

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_1~AT+edad+sexo, data=base, family=binomial)
# # plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ sexo + edad
## Model 2: alcohol_1 ~ AT + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     249.42                     
## 2       910     249.28  1  0.14283   0.7055

Se cumplen los supuestos del modelo, pero AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_1~AP+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ sexo + edad
## Model 2: alcohol_1 ~ AP + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     249.42                     
## 2       910     249.15  1   0.2707   0.6029

Se cumplen los supuestos del modelo, pero AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_1~AC+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ sexo + edad
## Model 2: alcohol_1 ~ AC + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       911     249.42                        
## 2       910     238.68  1   10.739 0.001049 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_1 ~ AC + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -20.3584     3.0927  -6.583 4.62e-11 ***
## AC            0.9601     0.3144   3.054  0.00226 ** 
## edad          0.7280     0.1404   5.185 2.16e-07 ***
## sexom         0.5152     0.3760   1.370  0.17063    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 238.68  on 910  degrees of freedom
## AIC: 246.68
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

Se cumplen los supuestos del modelo y AC fue significativo en la predicción de consumo de alcohol (\(OR\) = 2.61; \(IC_{95}\) = [1.45; 4.99]). De acuerdo a este resultado por cada unidad que aumenta la apertura al cambio, el riesgo de presentar una frecuencia de consumo elevada aumenta en un 161.19%.

Pred.: Conservación (C)

ajuste<-glm(alcohol_1~C+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad
## Model 2: alcohol_1 ~ C + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       912     250.18                        
## 2       910     239.05  2   11.127 0.003836 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_1 ~ C + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -12.4562     2.4353  -5.115 3.14e-07 ***
## C            -0.7435     0.2402  -3.095  0.00197 ** 
## edad          0.7153     0.1401   5.107 3.28e-07 ***
## sexom         0.3762     0.3743   1.005  0.31491    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 239.05  on 910  degrees of freedom
## AIC: 247.05
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

Se cumplen los supuestos del modelo y C fue significativo en la predicción de consumo de alcohol (\(OR\) = 0.48; \(IC_{95}\) = [0.29; 0.75]). De acuerdo a este resultado por cada unidad que aumenta la conservación, el riesgo de presentar una frecuencia de consumo elevada disminuye en un 52.45%.

Modelo final (Audit 1)

Se tomó como punto de base el modelo con edad, sexo y apertura al cambio por ser el que obtuvo una menor deviance y se comparó con el modelo que además tenía conservación.

mod_AC<-glm(alcohol_1~edad+sexo+AC, data=base, family=binomial)
mod_total <- glm(alcohol_1~edad+sexo+C+AC, data=base, family=binomial)
# plot(simulateResiduals(mod_total))
anova(mod_AC, mod_total, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_1 ~ edad + sexo + AC
## Model 2: alcohol_1 ~ edad + sexo + C + AC
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       910     238.68                          
## 2       909     222.25  1   16.432 5.042e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod_total)

## 
## Call:
## glm(formula = alcohol_1 ~ edad + sexo + C + AC, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -18.8828     3.2299  -5.846 5.03e-09 ***
## edad          0.7460     0.1481   5.038 4.70e-07 ***
## sexom         0.7513     0.3961   1.897 0.057855 .  
## C            -0.9223     0.2407  -3.831 0.000128 ***
## AC            1.2608     0.3378   3.732 0.000190 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.00  on 913  degrees of freedom
## Residual deviance: 222.25  on 909  degrees of freedom
## AIC: 232.25
## 
## Number of Fisher Scoring iterations: 7

OR_edad <- exp(mod_total$coefficients[2])
OR_sexo <- exp(mod_total$coefficients[3])
OR_C <- exp(mod_total$coefficients[4])
OR_AC <- exp(mod_total$coefficients[5])
IC_OR_edad<-exp(confint(mod_total))[2,]
IC_OR_sexo<-exp(confint(mod_total))[3,]
IC_OR_C<-exp(confint(mod_total))[4,]
IC_OR_AC<-exp(confint(mod_total))[5,]
R2 <- NagelkerkeR2(mod_total)$R2

Se cumplen los supuestos del modelo y edad, conservación y apertura al cambio fueron significativos en la predicción de consumo de alcohol (\(OR_{edad}\) = 2.11, \(IC_{95}\) = [1.6; 2.87]; \(OR_{Conservacion}\) = 0.4, \(IC_{95}\) = [0.24; 0.63]; \(OR_{AperturaAlCambio}\) = 3.53, \(IC_{95}\) = [1.88; 7.09]). El sexo obtuvo una significación marginal (\(OR_{sexo}\) = 2.12, \(IC_{95}\) = [0.98; 4.7]). De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar una frecuencia alta de consumo es 110.85% mayor por cada año que aumenta la edad, 252.84% mayor por cada unidad que aumenta la apertura al cambio, 111.97% mayor en varones que en mujeres, y 60.24% menor por cada unidad que aumenta la conservación.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.245) se estima que el modelo propuesto explica el 24.46% de la variabilidad de Audit 1.

Valores básicos

Dada la presencia de multicolinealidad evidenciada en la fase psicométrica, se decidió ajustar un modelo lineal generalizado (binomial) por cada valor por separado para indagar qué aspectos de la apertura al cambio y de la conservación son aquellos que mejor predicen la probabilidad de tener una frecuencia de consumo elevada.

mod_uni <- glm(alcohol_1~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_1~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_1~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_1~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_1~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_1~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_1~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_1~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_1~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_1~edad+sexo+seg, data=base, family=binomial)

# # plot(simulateResiduals(mod_seg))

tabla_valores <- data.frame(Valor=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"))
tabla_valores$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_valores$IC <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_valores$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_valores$R2 <- c(round(NagelkerkeR2(mod_uni)$R2,3),
                      round(NagelkerkeR2(mod_ben)$R2,3),
                      round(NagelkerkeR2(mod_adi)$R2,3),
                      round(NagelkerkeR2(mod_est)$R2,3),
                      round(NagelkerkeR2(mod_hed)$R2,3),
                      round(NagelkerkeR2(mod_log)$R2,3),
                      round(NagelkerkeR2(mod_pod)$R2,3),
                      round(NagelkerkeR2(mod_seg)$R2,3),
                      round(NagelkerkeR2(mod_con)$R2,3),
                      round(NagelkerkeR2(mod_tra)$R2,3))
tabla_valores$color <- ifelse(tabla_valores$p<.05&tabla_valores$OR<1,"#daaaaa",
                          ifelse(tabla_valores$p<.05&tabla_valores$OR>1,"#a5c3c6",NA))
tabla_valores$p <- ifelse(tabla_valores$p==0,"< .001", tabla_valores$p)



tabla <- kable(tabla_valores[1:5],
      "html",
      booktabs = T,
      align = c("c"),
      caption = "Predicción de probabilidad de frecuencia de consumo elevada según valores básicos") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_valores)){
  if(!is.na(tabla_valores$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_valores$color[i], color="#181414")
    
  }
}

tabla

Ninguno de los modelos presentó problemas de supuestos. Los tres valores de Apertura al cambio resultaron predictores significativos (controlando por sexo y edad), siendo la búsqueda de placer sensorial (hedonismo) y la búsqueda de independencia (autodirección) los OR más elevados (2.07 y 1.95 respectivamente). En cuanto a los valores de conservación, cumplir con las normas (conformidad) obtuvo el menor OR (0.55), seguido por la búsqueda de seguridad y estabilidad (seguridad; 0.68). El respeto y aceptación de las normas culturales (tradición) no resultó significativo.

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_2~AT+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AT + edad + sexo
##   Resid. Df Resid. Dev Df  Deviance Pr(>Chi)
## 1       911     431.34                      
## 2       910     431.34  1 2.488e-05    0.996

Se cumplen los supuestos del modelo, pero AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_2~AP+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AP + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     431.34                     
## 2       910     431.29  1 0.047318   0.8278

Se cumplen los supuestos del modelo, pero AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_2~AC+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ sexo + edad
## Model 2: alcohol_2 ~ AC + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     431.34                     
## 2       910     431.05  1  0.29079   0.5897

Se cumplen los supuestos del modelo, pero AC no resultó significativo en la predicción de consumo de alcohol.

Pred.: Conservación (C)

ajuste<-glm(alcohol_2~C+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_2 ~ edad
## Model 2: alcohol_2 ~ C + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       912     439.37                        
## 2       910     426.01  2   13.353  0.00126 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_2 ~ C + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -13.40995    1.72878  -7.757 8.70e-15 ***
## C            -0.35935    0.15800  -2.274  0.02295 *  
## edad          0.77434    0.09981   7.758 8.64e-15 ***
## sexom        -0.73349    0.27171  -2.700  0.00694 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 513.88  on 913  degrees of freedom
## Residual deviance: 426.01  on 910  degrees of freedom
## AIC: 434.01
## 
## Number of Fisher Scoring iterations: 6

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

Se cumplen los supuestos del modelo y C fue significativo en la predicción de presencia de CEEA (\(OR\) = 0.7; \(IC_{95}\) = [0.51; 0.95]). De acuerdo a este resultado por cada unidad que aumenta la conservación, el riesgo de presentar CEEA disminuye en un 30.19%.

Modelo final (Audit 2)

Se tomó como modelo final el que incluye sexo, edad y conservación por ser el único valor que mejoró significativamente las predicciones.

mod_total <- glm(alcohol_2~edad+sexo+C, data=base, family=binomial)
summary(mod_total)

## 
## Call:
## glm(formula = alcohol_2 ~ edad + sexo + C, family = binomial, 
##     data = base)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -13.40995    1.72878  -7.757 8.70e-15 ***
## edad          0.77434    0.09981   7.758 8.64e-15 ***
## sexom        -0.73349    0.27171  -2.700  0.00694 ** 
## C            -0.35935    0.15800  -2.274  0.02295 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 513.88  on 913  degrees of freedom
## Residual deviance: 426.01  on 910  degrees of freedom
## AIC: 434.01
## 
## Number of Fisher Scoring iterations: 6

OR_edad <- exp(mod_total$coefficients[2])
OR_sexo <- exp(mod_total$coefficients[3])
OR_C <- exp(mod_total$coefficients[4])
IC_OR_edad<-exp(confint(mod_total))[2,]
IC_OR_sexo<-exp(confint(mod_total))[3,]
IC_OR_C<-exp(confint(mod_total))[4,]
R2 <- NagelkerkeR2(mod_total)$R2

Se cumplen los supuestos del modelo y edad y conservación fueron significativos en la predicción de consumo de alcohol (\(OR_{edad}\) = 2.17, \(IC_{95}\) = [1.8; 2.66]; \(OR_{Conservacion}\) = 0.7, \(IC_{95}\) = [0.51; 0.95]). El sexo obtuvo no fue significativo. De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar CEEA es 116.92% mayor por cada año que aumenta la edad y 30.19% menor por cada unidad que aumenta la conservación.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.213) se estima que el modelo propuesto explica el 21.31% de la variabilidad de Audit 2.

Valores básicos

Dada la presencia de multicolinealidad evidenciada en la fase psicométrica, se decidió ajustar un modelo lineal generalizado (binomial) por cada valor por separado para indagar qué de la conservación son aquellos que mejor predicen la probabilidad de presentar CEEA.

mod_uni <- glm(alcohol_2~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_2~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_2~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_2~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_2~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_2~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_2~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_2~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_2~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_2~edad+sexo+seg, data=base, family=binomial)

# # plot(simulateResiduals(mod_con))

tabla_valores <- data.frame(Valor=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"))
tabla_valores$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_valores$IC <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_valores$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_valores$R2 <- c(round(NagelkerkeR2(mod_uni)$R2,3),
                      round(NagelkerkeR2(mod_ben)$R2,3),
                      round(NagelkerkeR2(mod_adi)$R2,3),
                      round(NagelkerkeR2(mod_est)$R2,3),
                      round(NagelkerkeR2(mod_hed)$R2,3),
                      round(NagelkerkeR2(mod_log)$R2,3),
                      round(NagelkerkeR2(mod_pod)$R2,3),
                      round(NagelkerkeR2(mod_seg)$R2,3),
                      round(NagelkerkeR2(mod_con)$R2,3),
                      round(NagelkerkeR2(mod_tra)$R2,3))
tabla_valores$color <- ifelse(tabla_valores$p<.05&tabla_valores$OR<1,"#daaaaa",
                          ifelse(tabla_valores$p<.05&tabla_valores$OR>1,"#a5c3c6",NA))
tabla_valores$p <- ifelse(tabla_valores$p==0,"< .001", tabla_valores$p)



tabla <- kable(tabla_valores[1:5],
      "html",
      booktabs = T,
      align = c("c"),
      caption = "Predicción de probabilidad de presencia de CEEA según valores básicos") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_valores)){
  if(!is.na(tabla_valores$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_valores$color[i], color="#181414")
    
  }
}

tabla

Ninguno de los modelos presentó problemas de supuestos. Los valores de conservación que resultaron significativos fueron seguridad (búsqueda de seguridad y estabilidad; OR = 0.76) y conformidad (cumplir con las normas y evitar disgustar a otrxs; OR = 0.78). El respeto y aceptación de las normas culturales (tradición) no resultó significativo.

Pred.: Autotrascendencia (AT)

ajuste<-glm(alcohol_3~AT+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ sexo + edad
## Model 2: alcohol_3 ~ AT + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     237.48                     
## 2       910     236.78  1  0.69984   0.4028

Se cumplen los supuestos del modelo, pero AT no resultó significativo en la predicción de consumo de alcohol.

Pred.: Autopromoción (AP)

ajuste<-glm(alcohol_3~AP+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ sexo + edad
## Model 2: alcohol_3 ~ AP + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       911     237.48                     
## 2       910     237.47  1 0.014582   0.9039

Se cumplen los supuestos del modelo, pero AP no resultó significativo en la predicción de consumo de alcohol.

Pred.: Apertura al cambio (AC)

ajuste<-glm(alcohol_3~AC+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_base, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ sexo + edad
## Model 2: alcohol_3 ~ AC + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       911     237.48                       
## 2       910     231.93  1   5.5514  0.01847 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_3 ~ AC + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.4544     3.3239  -6.755 1.42e-11 ***
## AC            0.6798     0.3033   2.241    0.025 *  
## edad          0.9639     0.1602   6.017 1.77e-09 ***
## sexom        -0.4917     0.3847  -1.278    0.201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 231.93  on 910  degrees of freedom
## AIC: 239.93
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

Se cumplen los supuestos del modelo y AC fue significativo en la predicción de consumo de alcohol (\(OR\) = 1.97; \(IC_{95}\) = [1.12; 3.68]). De acuerdo a este resultado por cada unidad que aumenta la apertura al cambio, el riesgo de presentar una frecuencia de CEEA elevada aumenta en un 97.34%.

Pred.: Conservación (C)

ajuste<-glm(alcohol_3~C+edad+sexo, data=base, family=binomial)
# plot(simulateResiduals(ajuste))
anova(mod_edad, ajuste, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad
## Model 2: alcohol_3 ~ C + edad + sexo
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       912     240.35                       
## 2       910     233.88  2   6.4744  0.03927 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(ajuste)

## 
## Call:
## glm(formula = alcohol_3 ~ C + edad + sexo, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -17.2944     2.8034  -6.169 6.87e-10 ***
## C            -0.4251     0.2279  -1.866   0.0621 .  
## edad          0.9593     0.1602   5.989 2.11e-09 ***
## sexom        -0.6245     0.3819  -1.635   0.1020    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 233.88  on 910  degrees of freedom
## AIC: 241.88
## 
## Number of Fisher Scoring iterations: 7

OR <- exp(ajuste$coefficients[2])
IC_OR<-exp(confint(ajuste))[2,]

Se cumplen los supuestos del modelo y C tuvo una significación marginal en la predicción de consumo de alcohol (\(OR\) = 0.65; \(IC_{95}\) = [0.41; 1.01]). De acuerdo a este resultado por cada unidad que aumenta la conservación, el riesgo de presentar una frecuencia de CEEA elevada disminuye en un 34.63%.

Modelo final (Audit 3)

Se tomó como punto de base el modelo con edad, sexo y apertura al cambio por ser el que obtuvo una menor deviance y se evaluó si la inclusión de conservación mejora las predicciones del modelo.

mod_AC<-glm(alcohol_3~edad+sexo+AC, data=base, family=binomial)
mod_total <- glm(alcohol_3~edad+sexo+C+AC, data=base, family=binomial)
# plot(simulateResiduals(mod_total))
anova(mod_AC, mod_total, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: alcohol_3 ~ edad + sexo + AC
## Model 2: alcohol_3 ~ edad + sexo + C + AC
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       910     231.93                       
## 2       909     226.28  1   5.6455   0.0175 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod_total)

## 
## Call:
## glm(formula = alcohol_3 ~ edad + sexo + C + AC, family = binomial, 
##     data = base)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -21.6106     3.4049  -6.347 2.20e-10 ***
## edad          0.9886     0.1651   5.988 2.13e-09 ***
## sexom        -0.4170     0.3919  -1.064  0.28725    
## C            -0.5278     0.2267  -2.328  0.01989 *  
## AC            0.8200     0.3160   2.595  0.00946 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 290.54  on 913  degrees of freedom
## Residual deviance: 226.28  on 909  degrees of freedom
## AIC: 236.28
## 
## Number of Fisher Scoring iterations: 7

OR_edad <- exp(mod_total$coefficients[2])
OR_sexo <- exp(mod_total$coefficients[3])
OR_C <- exp(mod_total$coefficients[4])
OR_AC <- exp(mod_total$coefficients[5])
IC_OR_edad<-exp(confint(mod_total))[2,]
IC_OR_sexo<-exp(confint(mod_total))[3,]
IC_OR_C<-exp(confint(mod_total))[4,]
IC_OR_AC<-exp(confint(mod_total))[5,]
R2 <- NagelkerkeR2(mod_total)$R2

Se cumplen los supuestos del modelo y la inclusión de Conservación mejora el ajuste significativamente. La edad, la conservación y la apertura al cambio fueron predictores significativos de la frecuencia alta de CEEA (\(OR_{edad}\) = 2.69, \(IC_{95}\) = [1.99; 3.8]; \(OR_{Conservacion}\) = 0.59, \(IC_{95}\) = [0.37; 0.91]; \(OR_{AperturaAlCambio}\) = 2.27, \(IC_{95}\) = [1.26; 4.35]). El sexo no resultó significativo. De acuerdo a este resultado, manteniendo fijas las demás variables, el riesgo de presentar una frecuencia alta de CEEA es 168.75% mayor por cada año que aumenta la edad, 127.05% mayor por cada unidad que aumenta la apertura al cambio, y 41.01% menor por cada unidad que aumenta la conservación.

A partir del pseudo-R2 de Nagelkerke (\(R^2\) = 0.249) se estima que el modelo propuesto explica el 24.93% de la variabilidad de Audit 3.

Valores básicos

Dada la presencia de multicolinealidad evidenciada en la fase psicométrica, se decidió ajustar un modelo lineal generalizado (binomial) por cada valor por separado para indagar qué aspectos de la apertura al cambio y de la conservación son aquellos que mejor predicen la probabilidad de tener una frecuencia de consumo elevada.

mod_uni <- glm(alcohol_3~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_3~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_3~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_3~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_3~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_3~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_3~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_3~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_3~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_3~edad+sexo+seg, data=base, family=binomial)

# # plot(simulateResiduals(mod_est))

tabla_valores <- data.frame(Valor=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"))
tabla_valores$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_valores$IC <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_valores$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_valores$R2 <- c(round(NagelkerkeR2(mod_uni)$R2,3),
                      round(NagelkerkeR2(mod_ben)$R2,3),
                      round(NagelkerkeR2(mod_adi)$R2,3),
                      round(NagelkerkeR2(mod_est)$R2,3),
                      round(NagelkerkeR2(mod_hed)$R2,3),
                      round(NagelkerkeR2(mod_log)$R2,3),
                      round(NagelkerkeR2(mod_pod)$R2,3),
                      round(NagelkerkeR2(mod_seg)$R2,3),
                      round(NagelkerkeR2(mod_con)$R2,3),
                      round(NagelkerkeR2(mod_tra)$R2,3))
tabla_valores$color <- ifelse(tabla_valores$p<.05&tabla_valores$OR<1,"#daaaaa",
                          ifelse(tabla_valores$p<.05&tabla_valores$OR>1,"#a5c3c6",NA))
tabla_valores$p <- ifelse(tabla_valores$p==0,"< .001", tabla_valores$p)



tabla <- kable(tabla_valores[1:5],
      "html",
      booktabs = T,
      align = c("c"),
      caption = "Predicción de probabilidad de frecuencia de CEEA elevada según valores básicos") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_valores)){
  if(!is.na(tabla_valores$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_valores$color[i], color="#181414")
    
  }
}

tabla

Ninguno de los modelos presentó problemas de supuestos. El único valor cuyo coeficiente fue significativo controlando por sexo y edad fue la búsqueda de placer sensorial (hedonismo; OR = 1.96). La búsqueda de seguridad y estabilidad (seguridad; OR = 0.76) y de novedad y desafíos (estimulación; OR = 1.43) obtuvieron una significación marginal en la predicción de la frecuencia de CEEA.

Edad

Como era de esperarse, en todos los modelos para predecir frecuencia de consumo, presencia de CEEA y frecuencia de CEEA la edad fue significativa, con ORs en los modelos finales de 2.11 para la frecuencia de consumo, 2.04 para la presencia de CEEA y 2.69 para la frecuencia de CEEA. De esta forma, por cada año que aumenta la edad, el riesgo de presentar una frecuencia de consumo elevada, de presencia de CEEA y de frecuencia de CEEA elevada aumentan más del doble (110.85%, 104.48% y 168.75% respectivamente).

Sexo

El sexo fue significativo sólo en el modelo de predicción de la presencia de CEEA, con un OR de 0.48 para los varones. Es decir, las chances de presentar un patrón de consumo de riesgo (CEEA) en varones es un 51.98% menor que en mujeres. En la última EMSE algo de esto aparece, aunque no recuerdo si específicamente con la presencia de CEEA.

En lo que respecta a la frecuencia de consumo, el sexo sólo alcanzó una significación marginal (p = .057) y con un OR en la dirección contraria (OR = 2.12). Pareciera ser que los varones tienen un riesgo 111.98% mayor que las mujeres de consumir alcohol con una frecuencia elevada, pero en lo que respecta a patrones de consumo ellos presentan un riesgo 51.98% menor que las adolescentes.

Valores

Los valores de autotrascendencia y de autopromoción no fueron predictores significativos en ninguno de los modelos.

Los valores de apertura al cambio permitieron predecir la probabilidad de presentar frecuencia alta de consumo (OR = 3.53) y frecuencia alta de CEEA (OR = 2.27), pero no presencia de CEEA. Más específicamente, que los motivos por los que los adolescentes se comportan como se comportan estén vinculados a la búsqueda de autonomía e independencia (autodirección; OR = 7.03), de novedad y desafíos (estimulación; OR = 4.81), y de disfrute y placer (hedonismo; 7.92) aumenta el riesgo de presentar una frecuencia de consumo. Por su parte, el riesgo de presentar una frecuencia de CEEA elevada aumenta un 609.93% por cada unidad que aumenta hedonismo (OR = 7.1), el único valor básico que predijo la frecuencia de CEEA.

Los valores de conservación predijeron significativamente la probabilidad de frecuencia de consumo (OR = 0.4), la presencia de CEEA (OR = 0.7) y la frecuencia de CEEA (OR = 0.59). Es decir, al tener niveles más altos de conservación, el riesgo de tener una frecuencia de consumo elevada, de presentar CEEA y de tener una frecuencia de CEEA elevada disminuye. En lo que respecta a los valores básicos, la búsqueda de seguridad y estabilidad (seguridad; OR = 0.51) y la evitación de incomodar a otros y la adecuación a las normas (conformidad; OR = 0.58) resultaron predictores significativos de la probabilidad de presentar una frecuencia de consumo elevada. Lo mismo ocurre en la predicción de la presencia de CEEA (\(OR_{seg}\) = 0.47; \(OR_{con}\) = 0.46). Sin embargo, ninguno de los modelos para predecir la frecuencia de CEEA con los valores básicos de conservación resultó significativo. Pareciera ser que la dimensión en su conjunto puede predecir el riesgo de presentar una frecuencia elevada de CEEA, pero no cada valor por separado.

Frecuencia de consumo

mod_base <- glm(alcohol_1~1, data=base, family=binomial)
mod <- glm(alcohol_1~edad+sexo+C+AC, data=base, family=binomial)

tabla_mod1 <- data.frame(Variable=c("Intercept","Age","Sex (male)","Conservation","Openness to Change"),
                         Beta=rep(NA,5),
                         SE=rep(NA,5),
                         p=rep(NA,5),
                         OR=rep(NA,5),
                         IC95=rep(NA,5),
                         Deviance=rep(NA,5),
                         df=rep(NA,5),
                         `p (model)`=rep(NA,5),
                         `Pseudo R2`=rep(NA,5),
                         check.names = F)
tabla_mod1$Beta <- round(mod$coefficients,2)
tabla_mod1$SE <- round(summary(mod)$coef[,2],2)
tabla_mod1$p <- round(summary(mod)$coef[,4],3)
tabla_mod1$OR <- round(exp(tabla_mod1$Beta),2)
tabla_mod1$IC95 <- paste0(round(exp(confint(mod))[,1],2)," - ",round(exp(confint(mod))[,2],2))
tabla_mod1$Deviance[1] <- round(anova(mod_base,mod, test = "Chisq")$`Deviance`[2],3)
tabla_mod1$df[1] <- anova(mod_base,mod, test = "Chisq")$`Df`[2]
tabla_mod1$`p (model)`[1] <- round(anova(mod_base,mod, test = "Chisq")$`Pr(>Chi)`[2],3)
tabla_mod1$`Pseudo R2` <- c("Cox Snell",
                            round(DescTools::PseudoR2(mod, c("CoxSnell")),3),
                            NA,
                            "Nagelkerke",
                            round(DescTools::PseudoR2(mod, c("Nagelkerke")),3))
tabla_mod1$p <- ifelse(tabla_mod1$p==0, "< .001", tabla_mod1$p)
tabla_mod1$`p (model)` <- ifelse(tabla_mod1$`p (model)`==0, "< .001", tabla_mod1$`p (model)`)


tabla <- kable(tabla_mod1,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de frecuencia de consumo elevada (Audit 1)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)


tabla

mod_base <- glm(alcohol_1~1, data=base, family=binomial)
mod_uni <- glm(alcohol_1~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_1~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_1~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_1~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_1~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_1~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_1~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_1~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_1~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_1~edad+sexo+seg, data=base, family=binomial)

tabla_mod1_10 <- data.frame(Variable=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"),
                         Beta=rep(NA,10),
                         SE=rep(NA,10),
                         p=rep(NA,10),
                         OR=rep(NA,10),
                         IC95=rep(NA,10),
                         Deviance=rep(NA,10),
                         df=rep(NA,10),
                         `p (model)`=rep(NA,10),
                         `Cox Snell R2`=rep(NA,10),
                         `Nagelkerke R2`=rep(NA,10),
                         check.names = F)
tabla_mod1_10$Beta <- c(round(mod_uni$coefficients[4],2),
                      round(mod_ben$coefficients[4],2),
                      round(mod_adi$coefficients[4],2),
                      round(mod_est$coefficients[4],2),
                      round(mod_hed$coefficients[4],2),
                      round(mod_log$coefficients[4],2),
                      round(mod_pod$coefficients[4],2),
                      round(mod_seg$coefficients[4],2),
                      round(mod_con$coefficients[4],2),
                      round(mod_tra$coefficients[4],2))
tabla_mod1_10$SE <- c(round(summary(mod_uni)$coef[4,2],2),
                      round(summary(mod_ben)$coef[4,2],2),
                      round(summary(mod_adi)$coef[4,2],2),
                      round(summary(mod_est)$coef[4,2],2),
                      round(summary(mod_hed)$coef[4,2],2),
                      round(summary(mod_log)$coef[4,2],2),
                      round(summary(mod_pod)$coef[4,2],2),
                      round(summary(mod_seg)$coef[4,2],2),
                      round(summary(mod_con)$coef[4,2],2),
                      round(summary(mod_tra)$coef[4,2],2))
tabla_mod1_10$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_mod1_10$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_mod1_10$IC95 <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_mod1_10$Deviance <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Deviance`[2],3))
tabla_mod1_10$df <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Df`[2],3))
tabla_mod1_10$`p (model)` <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Pr(>Chi)`[2],3))
tabla_mod1_10$`Cox Snell R2` <- c(round(DescTools::PseudoR2(mod_uni, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_est, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_log, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_con, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("CoxSnell")),3))
tabla_mod1_10$`Nagelkerke R2` <- c(round(DescTools::PseudoR2(mod_uni, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_est, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_log, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_con, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("Nagelkerke")),3))
tabla_mod1_10$color <- ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR<1,"#daaaaa",
                          ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR>1,"#a5c3c6",NA))
tabla_mod1_10$p <- ifelse(tabla_mod1_10$p==0, "< .001", tabla_mod1_10$p)
tabla_mod1_10$`p (model)` <- ifelse(tabla_mod1_10$`p (model)`==0, "< .001", tabla_mod1_10$`p (model)`)







tabla <- kable(tabla_mod1_10[1:10],
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de frecuencia de consumo elevada (Audit 1) según valores básicos (controlando por sexo y edad)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_mod1_10)){
  if(!is.na(tabla_mod1_10$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_mod1_10$color[i], color="#181414")
    
  }
}


tabla

Presencia de CEEA

mod_base <- glm(alcohol_2~1, data=base, family=binomial)
mod <- glm(alcohol_2~edad+sexo+C, data=base, family=binomial)

tabla_mod1 <- data.frame(Variable=c("Intercept","Age","Sex (male)","Conservation"),
                         Beta=rep(NA,4),
                         SE=rep(NA,4),
                         p=rep(NA,4),
                         OR=rep(NA,4),
                         IC95=rep(NA,4),
                         Deviance=rep(NA,4),
                         df=rep(NA,4),
                         `p (model)`=rep(NA,4),
                         `Pseudo R2`=rep(NA,4),
                         check.names = F)
tabla_mod1$Beta <- round(mod$coefficients,2)
tabla_mod1$SE <- round(summary(mod)$coef[,2],2)
tabla_mod1$p <- round(summary(mod)$coef[,4],3)
tabla_mod1$OR <- round(exp(tabla_mod1$Beta),2)
tabla_mod1$IC95 <- paste0(round(exp(confint(mod))[,1],2)," - ",round(exp(confint(mod))[,2],2))
tabla_mod1$Deviance[1] <- round(anova(mod_base,mod, test = "Chisq")$`Deviance`[2],3)
tabla_mod1$df[1] <- anova(mod_base,mod, test = "Chisq")$`Df`[2]
tabla_mod1$`p (model)`[1] <- round(anova(mod_base,mod, test = "Chisq")$`Pr(>Chi)`[2],3)
tabla_mod1$`Pseudo R2` <- c("Cox Snell",
                            round(DescTools::PseudoR2(mod, c("CoxSnell")),3),
                            "Nagelkerke",
                            round(DescTools::PseudoR2(mod, c("Nagelkerke")),3))
tabla_mod1$p <- ifelse(tabla_mod1$p==0, "< .001", tabla_mod1$p)
tabla_mod1$`p (model)` <- ifelse(tabla_mod1$`p (model)`==0, "< .001", tabla_mod1$`p (model)`)


tabla <- kable(tabla_mod1,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de presencia de CEEA (Audit 2)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)


tabla

mod_base <- glm(alcohol_2~1, data=base, family=binomial)
mod_uni <- glm(alcohol_2~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_2~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_2~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_2~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_2~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_2~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_2~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_2~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_2~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_2~edad+sexo+seg, data=base, family=binomial)

tabla_mod1_10 <- data.frame(Variable=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"),
                         Beta=rep(NA,10),
                         SE=rep(NA,10),
                         p=rep(NA,10),
                         OR=rep(NA,10),
                         IC95=rep(NA,10),
                         Deviance=rep(NA,10),
                         df=rep(NA,10),
                         `p (model)`=rep(NA,10),
                         `Cox Snell R2`=rep(NA,10),
                         `Nagelkerke R2`=rep(NA,10),
                         check.names = F)
tabla_mod1_10$Beta <- c(round(mod_uni$coefficients[4],2),
                      round(mod_ben$coefficients[4],2),
                      round(mod_adi$coefficients[4],2),
                      round(mod_est$coefficients[4],2),
                      round(mod_hed$coefficients[4],2),
                      round(mod_log$coefficients[4],2),
                      round(mod_pod$coefficients[4],2),
                      round(mod_seg$coefficients[4],2),
                      round(mod_con$coefficients[4],2),
                      round(mod_tra$coefficients[4],2))
tabla_mod1_10$SE <- c(round(summary(mod_uni)$coef[4,2],2),
                      round(summary(mod_ben)$coef[4,2],2),
                      round(summary(mod_adi)$coef[4,2],2),
                      round(summary(mod_est)$coef[4,2],2),
                      round(summary(mod_hed)$coef[4,2],2),
                      round(summary(mod_log)$coef[4,2],2),
                      round(summary(mod_pod)$coef[4,2],2),
                      round(summary(mod_seg)$coef[4,2],2),
                      round(summary(mod_con)$coef[4,2],2),
                      round(summary(mod_tra)$coef[4,2],2))
tabla_mod1_10$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_mod1_10$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_mod1_10$IC95 <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_mod1_10$Deviance <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Deviance`[2],3))
tabla_mod1_10$df <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Df`[2],3))
tabla_mod1_10$`p (model)` <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Pr(>Chi)`[2],3))
tabla_mod1_10$`Cox Snell R2` <- c(round(DescTools::PseudoR2(mod_uni, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_est, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_log, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_con, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("CoxSnell")),3))
tabla_mod1_10$`Nagelkerke R2` <- c(round(DescTools::PseudoR2(mod_uni, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_est, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_log, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_con, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("Nagelkerke")),3))
tabla_mod1_10$color <- ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR<1,"#daaaaa",
                          ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR>1,"#a5c3c6",NA))
tabla_mod1_10$p <- ifelse(tabla_mod1_10$p==0, "< .001", tabla_mod1_10$p)
tabla_mod1_10$`p (model)` <- ifelse(tabla_mod1_10$`p (model)`==0, "< .001", tabla_mod1_10$`p (model)`)







tabla <- kable(tabla_mod1_10[1:10],
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de presencia de CEEA (Audit 2) según valores básicos (controlando por sexo y edad)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_mod1_10)){
  if(!is.na(tabla_mod1_10$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_mod1_10$color[i], color="#181414")
    
  }
}


tabla

Frecuencia de CEEA

mod_base <- glm(alcohol_3~1, data=base, family=binomial)
mod <- glm(alcohol_3~edad+sexo+C+AC, data=base, family=binomial)

tabla_mod1 <- data.frame(Variable=c("Intercept","Age","Sex (male)","Conservation","Openness to Change"),
                         Beta=rep(NA,5),
                         SE=rep(NA,5),
                         p=rep(NA,5),
                         OR=rep(NA,5),
                         IC95=rep(NA,5),
                         Deviance=rep(NA,5),
                         df=rep(NA,5),
                         `p (model)`=rep(NA,5),
                         `Pseudo R2`=rep(NA,5),
                         check.names = F)
tabla_mod1$Beta <- round(mod$coefficients,2)
tabla_mod1$SE <- round(summary(mod)$coef[,2],2)
tabla_mod1$p <- round(summary(mod)$coef[,4],3)
tabla_mod1$OR <- round(exp(tabla_mod1$Beta),2)
tabla_mod1$IC95 <- paste0(round(exp(confint(mod))[,1],2)," - ",round(exp(confint(mod))[,2],2))
tabla_mod1$Deviance[1] <- round(anova(mod_base,mod, test = "Chisq")$`Deviance`[2],3)
tabla_mod1$df[1] <- anova(mod_base,mod, test = "Chisq")$`Df`[2]
tabla_mod1$`p (model)`[1] <- round(anova(mod_base,mod, test = "Chisq")$`Pr(>Chi)`[2],3)
tabla_mod1$`Pseudo R2` <- c("Cox Snell",
                            round(DescTools::PseudoR2(mod, c("CoxSnell")),3),
                            NA,
                            "Nagelkerke",
                            round(DescTools::PseudoR2(mod, c("Nagelkerke")),3))
tabla_mod1$p <- ifelse(tabla_mod1$p==0, "< .001", tabla_mod1$p)
tabla_mod1$`p (model)` <- ifelse(tabla_mod1$`p (model)`==0, "< .001", tabla_mod1$`p (model)`)


tabla <- kable(tabla_mod1,
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de frecuencia de CEEA elevada (Audit 3)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)


tabla

mod_base <- glm(alcohol_3~1, data=base, family=binomial)
mod_uni <- glm(alcohol_3~edad+sexo+uni, data=base, family=binomial)
mod_ben <- glm(alcohol_3~edad+sexo+ben, data=base, family=binomial)
mod_adi <- glm(alcohol_3~edad+sexo+adi, data=base, family=binomial)
mod_hed <- glm(alcohol_3~edad+sexo+hed, data=base, family=binomial)
mod_est <- glm(alcohol_3~edad+sexo+est, data=base, family=binomial)
mod_log <- glm(alcohol_3~edad+sexo+log, data=base, family=binomial)
mod_pod <- glm(alcohol_3~edad+sexo+pod, data=base, family=binomial)
mod_con <- glm(alcohol_3~edad+sexo+con, data=base, family=binomial)
mod_tra <- glm(alcohol_3~edad+sexo+tra, data=base, family=binomial)
mod_seg <- glm(alcohol_3~edad+sexo+seg, data=base, family=binomial)

tabla_mod1_10 <- data.frame(Variable=c("Universalismo",
                                    "Benevolencia",
                                    "Autodirección",
                                    "Estimulación",
                                    "Hedonismo",
                                    "Logro",
                                    "Poder",
                                    "Seguridad",
                                    "Conformidad",
                                    "Tradición"),
                         Beta=rep(NA,10),
                         SE=rep(NA,10),
                         p=rep(NA,10),
                         OR=rep(NA,10),
                         IC95=rep(NA,10),
                         Deviance=rep(NA,10),
                         df=rep(NA,10),
                         `p (model)`=rep(NA,10),
                         `Cox Snell R2`=rep(NA,10),
                         `Nagelkerke R2`=rep(NA,10),
                         check.names = F)
tabla_mod1_10$Beta <- c(round(mod_uni$coefficients[4],2),
                      round(mod_ben$coefficients[4],2),
                      round(mod_adi$coefficients[4],2),
                      round(mod_est$coefficients[4],2),
                      round(mod_hed$coefficients[4],2),
                      round(mod_log$coefficients[4],2),
                      round(mod_pod$coefficients[4],2),
                      round(mod_seg$coefficients[4],2),
                      round(mod_con$coefficients[4],2),
                      round(mod_tra$coefficients[4],2))
tabla_mod1_10$SE <- c(round(summary(mod_uni)$coef[4,2],2),
                      round(summary(mod_ben)$coef[4,2],2),
                      round(summary(mod_adi)$coef[4,2],2),
                      round(summary(mod_est)$coef[4,2],2),
                      round(summary(mod_hed)$coef[4,2],2),
                      round(summary(mod_log)$coef[4,2],2),
                      round(summary(mod_pod)$coef[4,2],2),
                      round(summary(mod_seg)$coef[4,2],2),
                      round(summary(mod_con)$coef[4,2],2),
                      round(summary(mod_tra)$coef[4,2],2))
tabla_mod1_10$p <- c(round(summary(mod_uni)$coefficients[4,4],3),
                     round(summary(mod_ben)$coefficients[4,4],3),
                     round(summary(mod_adi)$coefficients[4,4],3),
                     round(summary(mod_est)$coefficients[4,4],3),
                     round(summary(mod_hed)$coefficients[4,4],3),
                     round(summary(mod_log)$coefficients[4,4],3),
                     round(summary(mod_pod)$coefficients[4,4],3),
                     round(summary(mod_seg)$coefficients[4,4],3),
                     round(summary(mod_con)$coefficients[4,4],3),
                     round(summary(mod_tra)$coefficients[4,4],3))
tabla_mod1_10$OR <- c(round(exp(mod_uni$coefficients[4]),2),
                      round(exp(mod_ben$coefficients[4]),2),
                      round(exp(mod_adi$coefficients[4]),2),
                      round(exp(mod_est$coefficients[4]),2),
                      round(exp(mod_hed$coefficients[4]),2),
                      round(exp(mod_log$coefficients[4]),2),
                      round(exp(mod_pod$coefficients[4]),2),
                      round(exp(mod_seg$coefficients[4]),2),
                      round(exp(mod_con$coefficients[4]),2),
                      round(exp(mod_tra$coefficients[4]),2))
tabla_mod1_10$IC95 <- c(paste0(round(exp(confint(mod_uni))[4,1],2)," - ",round(exp(confint(mod_uni))[4,2],2)),
                      paste0(round(exp(confint(mod_ben))[4,1],2)," - ",round(exp(confint(mod_ben))[4,2],2)),
                      paste0(round(exp(confint(mod_adi))[4,1],2)," - ",round(exp(confint(mod_adi))[4,2],2)),
                      paste0(round(exp(confint(mod_est))[4,1],2)," - ",round(exp(confint(mod_est))[4,2],2)),
                      paste0(round(exp(confint(mod_hed))[4,1],2)," - ",round(exp(confint(mod_hed))[4,2],2)),
                      paste0(round(exp(confint(mod_log))[4,1],2)," - ",round(exp(confint(mod_log))[4,2],2)),
                      paste0(round(exp(confint(mod_pod))[4,1],2)," - ",round(exp(confint(mod_pod))[4,2],2)),
                      paste0(round(exp(confint(mod_seg))[4,1],2)," - ",round(exp(confint(mod_seg))[4,2],2)),
                      paste0(round(exp(confint(mod_con))[4,1],2)," - ",round(exp(confint(mod_con))[4,2],2)),
                      paste0(round(exp(confint(mod_tra))[4,1],2)," - ",round(exp(confint(mod_tra))[4,2],2)))
tabla_mod1_10$Deviance <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Deviance`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Deviance`[2],3))
tabla_mod1_10$df <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Df`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Df`[2],3))
tabla_mod1_10$`p (model)` <- c(round(anova(mod_base,mod_uni, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_ben, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_adi, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_est, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_hed, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_log, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_pod, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_seg, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_con, test = "Chisq")$`Pr(>Chi)`[2],3),
                            round(anova(mod_base,mod_tra, test = "Chisq")$`Pr(>Chi)`[2],3))
tabla_mod1_10$`Cox Snell R2` <- c(round(DescTools::PseudoR2(mod_uni, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_est, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_log, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_con, c("CoxSnell")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("CoxSnell")),3))
tabla_mod1_10$`Nagelkerke R2` <- c(round(DescTools::PseudoR2(mod_uni, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_ben, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_adi, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_est, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_hed, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_log, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_pod, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_seg, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_con, c("Nagelkerke")),3),
                                  round(DescTools::PseudoR2(mod_tra, c("Nagelkerke")),3))
tabla_mod1_10$color <- ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR<1,"#daaaaa",
                          ifelse(tabla_mod1_10$p<.05&tabla_mod1_10$OR>1,"#a5c3c6",NA))
tabla_mod1_10$p <- ifelse(tabla_mod1_10$p==0, "< .001", tabla_mod1_10$p)
tabla_mod1_10$`p (model)` <- ifelse(tabla_mod1_10$`p (model)`==0, "< .001", tabla_mod1_10$`p (model)`)







tabla <- kable(tabla_mod1_10[1:10],
      "html",
      booktabs = T,
      align = c("r"),
      caption = "Predicción de probabilidad de frecuencia de CEEA elevada (Audit 3) según valores básicos (controlando por sexo y edad)") %>%
  kable_styling(full_width = F,
                position = "center", font_size = 12)

for(i in 1:nrow(tabla_mod1_10)){
  if(!is.na(tabla_mod1_10$color[i])){
    tabla <- row_spec(tabla,i,background = tabla_mod1_10$color[i], color="#181414")
    
  }
}


tabla