El objetivo del estudio era el de evaluar si existía asociación entre la densidad mineral ósea (dmo) y las demás covariables.

Pregunta 1: IMC

Los investigadores deseaban evaluar primero si la densidad mineral ósea estaba asociada con el índice de masa corporal.

Se definen como categorías del índice de masa corporal

Normal: 18 a 24.9
Sobrepeso: 25 a 29.9
Obeso: 30 o mayor

Se categoriza la variable en los tres grupos

## 
##    Normal Sobrepeso  Obesidad 
##        94       194       111

Se construyen los gráficos de caja de la DMO para las tres categorías

## $Normal
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.00   61.25   73.00   73.28   84.00  136.00 
## 
## $Sobrepeso
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.00   62.00   72.00   72.71   82.75  130.00 
## 
## $Obesidad
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   35.00   65.00   77.00   77.06   88.00  121.00

La revisión del gráfico sugiere que el promedio de la DMO es similar entre los grupos.

Se identifican algunos valores atípicos. Para realizar la comparación de los promedios se utiliza un ANOVA. Se inicia con la validación de los supuestos de normalidad y homocesdasticidad.

## 
##  Shapiro-Wilk normality test
## 
## data:  normal$dmo
## W = 0.97392, p-value = 0.05029

## 
##  Shapiro-Wilk normality test
## 
## data:  sobrepeso$dmo
## W = 0.98551, p-value = 0.03895

## 
##  Shapiro-Wilk normality test
## 
## data:  obesidad$dmo
## W = 0.98904, p-value = 0.5595

En el grupo de sobrepeso se rechaza la hipótesis nula de normalidad. Se continua con la valoración de la homocesdasticidad

## 
##  F test to compare two variances
## 
## data:  normal$dmo and sobrepeso$dmo
## F = 1.0642, num df = 96, denom df = 198, p-value = 0.7077
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7606379 1.5206858
## sample estimates:
## ratio of variances 
##           1.064233

## 
##  F test to compare two variances
## 
## data:  normal$dmo and obesidad$dmo
## F = 0.83892, num df = 96, denom df = 103, p-value = 0.3846
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5656741 1.2476178
## sample estimates:
## ratio of variances 
##          0.8389177

## 
##  F test to compare two variances
## 
## data:  sobrepeso$dmo and obesidad$dmo
## F = 0.78828, num df = 198, denom df = 103, p-value = 0.1567
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5569406 1.0955227
## sample estimates:
## ratio of variances 
##          0.7882841

En ninguna comparación se rechaza la hipótesis nula de homogeneidad de las varianzas. Se realiza entonces la comparación con el ANOVA

modelo<-lm(osteo$dmo~osteo$imccat, data=osteo)
anova(modelo)

## Analysis of Variance Table
## 
## Response: osteo$dmo
##               Df Sum Sq Mean Sq F value  Pr(>F)  
## osteo$imccat   2   1412  705.79  2.3186 0.09974 .
## Residuals    396 120543  304.40                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El valor F obtenido es de 2.3186 y no se rechaza la hipótesis nula, lo cual indica que no hay diferencia en el promedio entre los grupos.

Pregunta 2: Edad

Los investigadores también deseaban evaluar si la densidad mineral ósea estaba asociada con la edad.

Se categoriza la edad en tres grupos (<50, de 50-59 y >59 años)

## 
##     <50 50 a 59     >59 
##     173     157      70

## # A tibble: 3 x 2
##   catedad `mean(dmo)`
##   <fct>         <dbl>
## 1 <50            80.3
## 2 50 a 59        72.3
## 3 >59            62.8

Se construye un gráfico de caja con las categorías construídas para la variable DMO

boxplot(tapply(osteo$dmo, osteo$catedad, summary), main="DMO por grupo de edad")

El gráfico de caja sugiere que hay diferencias entre el promedio de densidad mineral ósea entre los grupos. Se realiza entonces una comparación mediante un ANOVA. Se inicia con la validación de los supuestos de normalidad y homocedasticidad.

library(tidyverse)
Menor50<-osteo %>% filter(edad<50)
De50a59<-osteo %>% filter(edad>=50 & edad<59)
Mayor59<-osteo %>% filter(edad>=59)

Validación del supuesto de normalidad de la variable DMO en cada grupo de edad

shapiro.test(Menor50$dmo)

## 
##  Shapiro-Wilk normality test
## 
## data:  Menor50$dmo
## W = 0.96882, p-value = 0.001742

shapiro.test(De50a59$dmo)

## 
##  Shapiro-Wilk normality test
## 
## data:  De50a59$dmo
## W = 0.98185, p-value = 0.03112

shapiro.test(Mayor59$dmo)

## 
##  Shapiro-Wilk normality test
## 
## data:  Mayor59$dmo
## W = 0.95272, p-value = 0.003039

La prueba de normalidad rechaza la hipótesis nula de normalidad en los tres grupos. Se valora mediante un histograma la distribución.

azult<-rgb(0,0,155, alpha = 30, maxColorValue = 155)
rojot<-rgb(155,0,0, alpha = 30, maxColorValue = 155)
verdet<-rgb(0,155,0, alpha = 30, maxColorValue = 155)

par(mfrow=c(1,3))
hist(Menor50$dmo, col=azult, main="DMO en < 50", xlab="DMO", ylab="Frecuencia")
hist(De50a59$dmo, col=rojot, main="DMO en 50 a 59",xlab="DMO", ylab="Frecuencia")
hist(Mayor59$dmo, col=verdet, main="DMO en > 59",xlab="DMO", ylab="Frecuencia")

Validación del supuesto de homocedasticidad

var.test(Menor50$dmo, De50a59$dmo)

## 
##  F test to compare two variances
## 
## data:  Menor50$dmo and De50a59$dmo
## F = 1.0236, num df = 149, denom df = 162, p-value = 0.8827
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7472418 1.4052322
## sample estimates:
## ratio of variances 
##           1.023638

var.test(Menor50$dmo, Mayor59$dmo)

## 
##  F test to compare two variances
## 
## data:  Menor50$dmo and Mayor59$dmo
## F = 0.99129, num df = 149, denom df = 86, p-value = 0.95
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6731121 1.4319442
## sample estimates:
## ratio of variances 
##          0.9912922

var.test(De50a59$dmo, Mayor59$dmo)

## 
##  F test to compare two variances
## 
## data:  De50a59$dmo and Mayor59$dmo
## F = 0.9684, num df = 162, denom df = 86, p-value = 0.8498
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6605728 1.3895803
## sample estimates:
## ratio of variances 
##          0.9684013

No se rechaza la hipótesis nula de homogeneidad de las varianzas. Se procede entonces con la construcción del ANOVA

modelo2<-lm(osteo$dmo~osteo$catedad, data=osteo)
anova(modelo2)

## Analysis of Variance Table
## 
## Response: osteo$dmo
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## osteo$catedad   2  16000  8000.2  29.904 7.991e-13 ***
## Residuals     397 106208   267.5                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El valor del estadístico F es de 29.904 que corresponde a un valor menor de 0.00001. Se rechaza la hipótesis nula de que los promedios son iguales a favor de la hipótesis alternativa de que al menos uno es diferente. Para identificar entre grupos se encuentran las diferencias, se realizar una prueba de Tukey.

modelo_tukey<-aov(osteo$dmo~osteo$catedad, data = osteo)
TukeyHSD(modelo_tukey)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = osteo$dmo ~ osteo$catedad, data = osteo)
## 
## $`osteo$catedad`
##                   diff       lwr        upr     p adj
## 50 a 59-<50  -7.971724 -12.21307  -3.730377 0.0000376
## >59-<50     -17.448885 -22.89957 -11.998204 0.0000000
## >59-50 a 59  -9.477161 -15.00727  -3.947051 0.0001959

Se encuentra que el promedio es diferente entre todos los grupos comparados con valores significativos menores a 0,001 en todos los casos.

Pregunta 3: Correlación

¿Existe correlación entre la variable desenlace (dmo) y las variables independientes? Se construye una matriz de correlación entre las variables cuantitativas.

Se construye una matriz con los valores p de las correlaciones

Se presenta gráficamente la información de la matríz de correlación

Se encuentra que la variable DMO tiene correlación negativa de -0.37 con la edad y en menor grado (-0.13) con la edad de la menarquia. La correlación entre la DMO y la edad es significativa estadísticamente con un valor p de 0.02.

Punto 4: Modelo de asociación

Para construir el modelo se inicia con un modelo lineal simple entre la variable dmo y la variable edad que fue la que tuvo la correlación mas alta. Se declaran además como factores las variables categóricas antes de su inclusión en el modelo.

Se seleccionan los casos con datos completos

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	125.2	7.373	16.98	7.788e-47
osteoCom$edad	-0.9576	0.1391	-6.883	2.986e-11

Fitting linear model: osteoCom$dmo ~ osteoCom$edad
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
330	16.58	0.1262	0.1235

Analysis of Variance Table
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
osteoCom$edad	1	13026	13026	47.38	2.986e-11
Residuals	328	90182	274.9	NA	NA

## Analysis of Variance Table
## 
## Response: osteoCom$dmo
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## osteoCom$edad   1  13026 13026.0  47.377 2.986e-11 ***
## Residuals     328  90182   274.9                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En el modelo linear entre la DMO y la edad, se encuentra asociación. El modelo se interpreta como que, por cada año de incremento de la edad, se reduce la densidad mineral ósea en 0.9576 unidades. De los casos incluidos, fueron seleccionados para el modelo 330
Se construye un modelo completo con todas las variables para evaluar su comportamiento.

library(pander)
pander(anova(lm(osteoCom$dmo~osteoCom$edad+osteoCom$peso+osteoCom$talla+
           osteoCom$menarquia+osteoCom$menopausia+osteoCom$antmaternos+
           osteoCom$activ_laboral+osteoCom$ciclomenstrual+osteoCom$nivel_ed+
           osteoCom$paridad+osteoCom$numerohijos+osteoCom$abortos+
           osteoCom$imccat)))

Analysis of Variance Table
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
osteoCom$edad	1	13026	13026	48.58	1.905e-11
osteoCom$peso	1	3326	3326	12.4	0.0004929
osteoCom$talla	1	4.986	4.986	0.01859	0.8916
osteoCom$menarquia	1	799.6	799.6	2.982	0.08519
osteoCom$menopausia	2	55.38	27.69	0.1033	0.9019
osteoCom$antmaternos	1	38.13	38.13	0.1422	0.7064
osteoCom$activ_laboral	1	226.1	226.1	0.8432	0.3592
osteoCom$ciclomenstrual	1	3.4	3.4	0.01268	0.9104
osteoCom$nivel_ed	4	986	246.5	0.9193	0.4529
osteoCom$paridad	1	490	490	1.827	0.1774
osteoCom$numerohijos	1	14.94	14.94	0.05573	0.8135
osteoCom$abortos	1	256.1	256.1	0.955	0.3292
osteoCom$imccat	2	588	294	1.096	0.3354
Residuals	311	83393	268.1	NA	NA

La variables edad y peso son estadísticamente significativas en este modelo. La variable menarquia estuvo cerca al límite de la significancia. Las demás variables no reportan valores de p significativos.
Se decide construir un modelo ordenado, iniciando con las variables que tuvieron la correlación más alta. El modelo inicial incluyó solamente la edad como variable predictora.
Se construyen las correlaciones parciales, iniciando por la inclusión de la variable menarquia que fue la siguiente con mayor importancia en la matriz de correlación.

library(ppcor)
# Primero la que se quiere evaluar, luego la independiente, luego la incluida
pander(pcor.test(osteoCom$menarquia, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
-0.1085	0.04932	-1.973	330	1	pearson

pander(pcor.test(osteoCom$peso, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.192	0.0004607	3.538	330	1	pearson

pander(pcor.test(osteoCom$imccat, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.1795	0.001073	3.3	330	1	pearson

pander(pcor.test(osteoCom$menopausia, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.01478	0.7894	0.2673	330	1	pearson

pander(pcor.test(osteoCom$nivel_ed, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
-0.08808	0.1108	-1.599	330	1	pearson

pander(pcor.test(osteoCom$paridad, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.09578	0.0828	1.74	330	1	pearson

pander(pcor.test(osteoCom$antmaternos, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.01847	0.7386	0.334	330	1	pearson

pander(pcor.test(osteoCom$numerohijos, osteoCom$dmo,osteoCom[,c("edad")] ))

estimate	p.value	statistic	n	gp	Method
0.1001	0.06974	1.82	330	1	pearson

La segunda variable a incluir en el modelo es el peso, el cual tuvo el valor de p mas bajo.

modelo_tercero<-lm(osteoCom$dmo~osteoCom$edad+osteoCom$peso, data=osteoCom)
pander(summary(modelo_tercero))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	108.9	8.575	12.7	2.483e-30
osteoCom$edad	-1.032	0.1383	-7.458	7.99e-13
osteoCom$peso	0.2949	0.08334	3.538	0.0004607

Fitting linear model: osteoCom$dmo ~ osteoCom$edad + osteoCom$peso
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
330	16.3	0.1584	0.1533

pander(anova(modelo_tercero))

Analysis of Variance Table
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
osteoCom$edad	1	13026	13026	49.04	1.432e-11
osteoCom$peso	1	3326	3326	12.52	0.0004607
Residuals	327	86856	265.6	NA	NA

La inclusión de la nueva variable mejora el ajuste del modelo, con aumento del valor del $R^2$.

Se construyen nuevamente las correlaciones parciales. La primera variable que se evalua es el IMC categorizado, seguido de la menarquia por ser significativas en la correlación parcial anterior.

pander(pcor.test(osteoCom$imccat, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
0.05147	0.3527	0.9306	330	2	pearson

pander(pcor.test(osteoCom$menarquia, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
-0.09394	0.08941	-1.704	330	2	pearson

pander(pcor.test(osteoCom$numerohijos, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
0.05321	0.3367	0.9622	330	2	pearson

pander(pcor.test(osteoCom$paridad, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
0.08217	0.1376	1.489	330	2	pearson

pander(pcor.test(osteoCom$nivel_ed, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
-0.04596	0.4067	-0.8308	330	2	pearson

pander(pcor.test(osteoCom$antmaternos, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
0.01027	0.853	0.1855	330	2	pearson

pander(pcor.test(osteoCom$menopausia, osteoCom$dmo,osteoCom[,c("edad", "peso")] ))

estimate	p.value	statistic	n	gp	Method
0.01239	0.8231	0.2238	330	2	pearson

En la correlación parcial, el IMC categorizado no fue estadísticamente significativo. La menarquia estuvo en el límite de la significancia y se decide evaluar el efecto de su inclusión en el modelo.

modelo_cuarto<-lm(osteoCom$dmo~osteoCom$edad+osteoCom$peso+osteoCom$menarquia, data=osteoCom)
pander(summary(modelo_cuarto))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	121.7	11.37	10.7	4.154e-23
osteoCom$edad	-1.009	0.1386	-7.277	2.559e-12
osteoCom$peso	0.2826	0.08341	3.388	0.0007906
osteoCom$menarquia	-1.035	0.6074	-1.704	0.08941

Fitting linear model: osteoCom$dmo ~ osteoCom$edad + osteoCom$peso + osteoCom$menarquia
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
330	16.25	0.1659	0.1582

anova(modelo_cuarto)

## Analysis of Variance Table
## 
## Response: osteoCom$dmo
##                     Df Sum Sq Mean Sq F value    Pr(>F)    
## osteoCom$edad        1  13026 13026.0 49.3263 1.268e-11 ***
## osteoCom$peso        1   3326  3325.7 12.5937 0.0004439 ***
## osteoCom$menarquia   1    766   766.4  2.9022 0.0894130 .  
## Residuals          326  86090   264.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pander(anova(modelo_cuarto))

Analysis of Variance Table
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
osteoCom$edad	1	13026	13026	49.33	1.268e-11
osteoCom$peso	1	3326	3326	12.59	0.0004439
osteoCom$menarquia	1	766.4	766.4	2.902	0.08941
Residuals	326	86090	264.1	NA	NA

Con las tres variables incluidas el $R^2$ aumenta a 15.82. Se considera entonces que el modelo final es:

$dmo=\beta_0+\beta_1 edad+\beta_2peso+\beta_3 menarquia$

La tabla ANOVA final que se construye con los valores correctos es:

         Df Sum Sq Mean Sq F value    Pr(>F)   
   edad   1  13026 13026.0  47.377 2.986e-11 *** 
   peso   1   3326  3325.7  12.521 0.0004607 ***

menarquia 1 766 766.4 2.9022 0.0894130

Diagnóstico del modelo

Se realiza primero la identificación del mejor modelo para proceder a su evaluación. En el modelo propuesto fueron incluidas tres variables que fueron edad, peso y menarquia. Se construyen entonces dos modelos: uno incluyendo las tres variables y otro excluyendo la variable menarquia que no tuvo un valor p significativo en la correlación parcial cuando se encontraban incluidas las otras dos variables.

# El modelo con las tres variables se llama modelocuarto 
modelo_quinto<-lm(osteoCom$dmo~osteoCom$edad+osteoCom$peso)
anova(modelo_quinto)

## Analysis of Variance Table
## 
## Response: osteoCom$dmo
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## osteoCom$edad   1  13026 13026.0  49.041 1.432e-11 ***
## osteoCom$peso   1   3326  3325.7  12.521 0.0004607 ***
## Residuals     327  86856   265.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(modelo_cuarto)

## Analysis of Variance Table
## 
## Response: osteoCom$dmo
##                     Df Sum Sq Mean Sq F value    Pr(>F)    
## osteoCom$edad        1  13026 13026.0 49.3263 1.268e-11 ***
## osteoCom$peso        1   3326  3325.7 12.5937 0.0004439 ***
## osteoCom$menarquia   1    766   766.4  2.9022 0.0894130 .  
## Residuals          326  86090   264.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(modelo_cuarto, modelo_quinto)

## Analysis of Variance Table
## 
## Model 1: osteoCom$dmo ~ osteoCom$edad + osteoCom$peso + osteoCom$menarquia
## Model 2: osteoCom$dmo ~ osteoCom$edad + osteoCom$peso
##   Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
## 1    326 86090                              
## 2    327 86856 -1    -766.4 2.9022 0.08941 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No se encuentran diferencias significativas entre los dos modelos, incluyendo o no a la menarquia por lo que se excluye. El modelo final sobre el cual se hace el diagnóstico incluye únicamente la variable edad y la variable peso. Para el diagnóstico del modelo se realizan los siguientes análisis: + Presencia de valores atípicos + Validación de los supuestos + Evaluación de colinealidad

Presencia de valores atípicos

Se inicia realizado el análisis de las residuales del modelo. Se utiliza inicialmente la valoración gráfica. La gráfica de los valores residuales crudas muestra una distribución simétrica alrededor del cero. El gráfico QQ normal de las residuales estandarizadas parece identificar como valores atípicos a tres valores: 53, 272 y 176. La evaluación gráfica de la distancia de Cook no parece identificar valores atípicos.

plot(modelo_quinto)

Se realizan pruebas para identificar valores atipicos asi

studentizadas<-rstudent(modelo_quinto)
studentizadas[which(abs(studentizadas)>2)]

##        25        53       104       152       176       194       197       240 
##  3.057240 -4.294231  2.018241 -2.696126  3.646492  2.441851 -2.318178  2.954056 
##       268       272       287       292       293       298 
##  2.540889  3.473689  2.426025  2.203394  2.328830  2.130944

Se identifican varios valores atípicos adicionalmente a lo que se había identificado gráficamente. Ahora se busca identificar la presencia de observaciones influyentes.

summary(influence.measures(modelo_quinto))

## Potentially influential observations of
##   lm(formula = osteoCom$dmo ~ osteoCom$edad + osteoCom$peso) :
## 
##     dfb.1_ dfb.ostCm$d dfb.ostCm$p dffit   cov.r   cook.d hat    
## 10  -0.02   0.03       -0.01        0.03    1.03_*  0.00   0.02  
## 16   0.02  -0.02        0.00       -0.03    1.03_*  0.00   0.02  
## 25  -0.17   0.12        0.12        0.25    0.93_*  0.02   0.01  
## 38   0.03   0.03       -0.08       -0.09    1.04_*  0.00   0.03_*
## 53  -0.17   0.16        0.01       -0.29_*  0.86_*  0.03   0.00  
## 65   0.05   0.04       -0.14       -0.15    1.08_*  0.01   0.07_*
## 81  -0.03   0.01        0.04        0.05    1.04_*  0.00   0.03_*
## 89   0.01   0.06       -0.09       -0.11    1.03_*  0.00   0.03  
## 97   0.04  -0.08        0.04       -0.09    1.03_*  0.00   0.02  
## 122  0.04   0.00       -0.06       -0.06    1.04_*  0.00   0.03_*
## 152  0.13   0.00       -0.23       -0.27    0.95_*  0.02   0.01  
## 171  0.07  -0.09        0.01       -0.10    1.03_*  0.00   0.02  
## 176  0.24  -0.15       -0.15        0.30_*  0.90_*  0.03   0.01  
## 181  0.29  -0.17       -0.25       -0.34_*  1.00    0.04   0.03_*
## 191  0.08   0.10       -0.26       -0.28    1.02    0.03   0.03_*
## 194 -0.01  -0.05        0.10        0.17    0.96_*  0.01   0.00  
## 196  0.00   0.01       -0.02       -0.02    1.04_*  0.00   0.03_*
## 216  0.03   0.04       -0.10       -0.10    1.04_*  0.00   0.03_*
## 236  0.04   0.03       -0.11       -0.11    1.04_*  0.00   0.03_*
## 240  0.11  -0.17        0.07        0.24    0.94_*  0.02   0.01  
## 253  0.04  -0.06        0.01       -0.07    1.03_*  0.00   0.02  
## 268  0.04  -0.06        0.03        0.15    0.96_*  0.01   0.00  
## 272 -0.49   0.41        0.24        0.55_*  0.93_*  0.10   0.02  
## 286  0.04  -0.03       -0.03       -0.05    1.04_*  0.00   0.03_*
## 287  0.19  -0.14       -0.09        0.22    0.96_*  0.02   0.01  
## 292  0.12  -0.14        0.01        0.19    0.97_*  0.01   0.01

Se identifican múltiples observaciones potencialmente influyentes.

library(car)

## Warning: package 'car' was built under R version 3.6.3

## Loading required package: carData

## 
## Attaching package: 'car'

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

## The following object is masked from 'package:purrr':
## 
##     some

outlierTest(modelo_quinto)

##     rstudent unadjusted p-value Bonferroni p
## 53 -4.294231         2.3155e-05    0.0076413

La prueba de bonferroni identifica como valor atípico el valor 53.

Validación de supuestos del modelo

Homocedasticidad

library(lmtest)

## Warning: package 'lmtest' was built under R version 3.6.3

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(modelo_quinto)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_quinto
## BP = 1.0357, df = 2, p-value = 0.5958

No se rechaza la hipotésis nula de que los datos son homocedasticos. Independencia

dwtest(modelo_quinto)

## 
##  Durbin-Watson test
## 
## data:  modelo_quinto
## DW = 2.0807, p-value = 0.7675
## alternative hypothesis: true autocorrelation is greater than 0

No se rechaza la hipótesis nula de que los errores no están correlacionados por lo que se valida la independencia Linealidad

crPlots(modelo_quinto)

La gráfica de edad muestra linealidad pero se aprecian desviaciones en la variable peso. Pruebas globales

library(gvlma)
gvmodel <- gvlma(modelo_quinto) 
summary(modelo_quinto)

## 
## Call:
## lm(formula = osteoCom$dmo ~ osteoCom$edad + osteoCom$peso)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -68.037 -10.730  -1.653  10.802  58.143 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   108.94716    8.57542  12.705  < 2e-16 ***
## osteoCom$edad  -1.03167    0.13834  -7.458 7.99e-13 ***
## osteoCom$peso   0.29489    0.08334   3.538 0.000461 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.3 on 327 degrees of freedom
## Multiple R-squared:  0.1584, Adjusted R-squared:  0.1533 
## F-statistic: 30.78 on 2 and 327 DF,  p-value: 5.647e-13

library(olsrr)

## Warning: package 'olsrr' was built under R version 3.6.3

## 
## Attaching package: 'olsrr'

## The following object is masked from 'package:MASS':
## 
##     cement

## The following object is masked from 'package:datasets':
## 
##     rivers

ols_plot_resid_stud_fit(modelo_quinto)

ols_plot_cooksd_bar(modelo_quinto)

El gráfico define una distancia de Cook de 0.012 como el limite para considerar un valor como atípico. Aunque varias observaciones están por encima del umbral, solamente una parece ser un valor influyente.

Colinearidad

summary(influence.measures(modelo_quinto))

## Potentially influential observations of
##   lm(formula = osteoCom$dmo ~ osteoCom$edad + osteoCom$peso) :
## 
##     dfb.1_ dfb.ostCm$d dfb.ostCm$p dffit   cov.r   cook.d hat    
## 10  -0.02   0.03       -0.01        0.03    1.03_*  0.00   0.02  
## 16   0.02  -0.02        0.00       -0.03    1.03_*  0.00   0.02  
## 25  -0.17   0.12        0.12        0.25    0.93_*  0.02   0.01  
## 38   0.03   0.03       -0.08       -0.09    1.04_*  0.00   0.03_*
## 53  -0.17   0.16        0.01       -0.29_*  0.86_*  0.03   0.00  
## 65   0.05   0.04       -0.14       -0.15    1.08_*  0.01   0.07_*
## 81  -0.03   0.01        0.04        0.05    1.04_*  0.00   0.03_*
## 89   0.01   0.06       -0.09       -0.11    1.03_*  0.00   0.03  
## 97   0.04  -0.08        0.04       -0.09    1.03_*  0.00   0.02  
## 122  0.04   0.00       -0.06       -0.06    1.04_*  0.00   0.03_*
## 152  0.13   0.00       -0.23       -0.27    0.95_*  0.02   0.01  
## 171  0.07  -0.09        0.01       -0.10    1.03_*  0.00   0.02  
## 176  0.24  -0.15       -0.15        0.30_*  0.90_*  0.03   0.01  
## 181  0.29  -0.17       -0.25       -0.34_*  1.00    0.04   0.03_*
## 191  0.08   0.10       -0.26       -0.28    1.02    0.03   0.03_*
## 194 -0.01  -0.05        0.10        0.17    0.96_*  0.01   0.00  
## 196  0.00   0.01       -0.02       -0.02    1.04_*  0.00   0.03_*
## 216  0.03   0.04       -0.10       -0.10    1.04_*  0.00   0.03_*
## 236  0.04   0.03       -0.11       -0.11    1.04_*  0.00   0.03_*
## 240  0.11  -0.17        0.07        0.24    0.94_*  0.02   0.01  
## 253  0.04  -0.06        0.01       -0.07    1.03_*  0.00   0.02  
## 268  0.04  -0.06        0.03        0.15    0.96_*  0.01   0.00  
## 272 -0.49   0.41        0.24        0.55_*  0.93_*  0.10   0.02  
## 286  0.04  -0.03       -0.03       -0.05    1.04_*  0.00   0.03_*
## 287  0.19  -0.14       -0.09        0.22    0.96_*  0.02   0.01  
## 292  0.12  -0.14        0.01        0.19    0.97_*  0.01   0.01

Se identifican varias observaciones potencialmente influyentes.

ols_vif_tol(modelo_quinto)

##       Variables Tolerance      VIF
## 1 osteoCom$edad 0.9771211 1.023415
## 2 osteoCom$peso 0.9771211 1.023415

El factor de inflación de la varianza es menor a 10. Se considera por este criterio que no hay multicolinealidad

ols_coll_diag(modelo_quinto)

## Tolerance and Variance Inflation Factor
## ---------------------------------------
##       Variables Tolerance      VIF
## 1 osteoCom$edad 0.9771211 1.023415
## 2 osteoCom$peso 0.9771211 1.023415
## 
## 
## Eigenvalue and Condition Index
## ------------------------------
##    Eigenvalue Condition Index   intercept osteoCom$edad osteoCom$peso
## 1 2.975109147         1.00000 0.001230207   0.001678083   0.002711659
## 2 0.017904185        12.89063 0.028140586   0.275647826   0.846825711
## 3 0.006986668        20.63557 0.970629207   0.722674090   0.150462631

Ninguno de los índices de condición es mayor a 30 por lo que se considera que no hay problemas de colinealidad.

Interpretación del diagnóstico

El modelo construido como un modelo de asociación tiene pocos valores atípicos y algunas obsrvaciones influyentes. Fue posible validar los supuestos del modelo y no se detecta colinealidad. La relación entre la variable edad y la DMO tiene un aspecto lineal pero la variable peso se desvia de este supuesto. Aunque pareciera razonable utilizar un modelo lineal, es posible que otro tipo de modelamiento pueda explicar mejor la relación entre las variables.

Taller1

Leonardo Enciso

28/4/2020