Clase II_Maestria

Chi-Cuadrado

El estadistico de chi-cuadrado, es util cuando estamos trabajando con variables nominales, categoricas o estamos trabajando con alguna clase de clasificacion Su formula estadistica esta dada por:

Donde: O= se refiere a las frecuencias observadas E= frecuencias esperadas

Requerimientos para trabajar con la distribucion de chi-cuadrado:

La muestra debe ser tomadas al azar; Variables medidas deben ser independientes; Los datos deben ser reportados en frecuencias absolutas (no porcentajes); Valores / categorias de variables dependientes e independientes deben ser mutuamente excluyentes;

El estadistico de chi resulta util para trabajar datos, como frecuencias, o datos categoricos.

Ejemplo 1. Hay 15 hombres y 19 mujeres en un aula de clase. Como ven en los datos, no nos interesa cada dato como independiente, si no como su frecuencia de la clase o categoria dada, en este caso la clase es el genero (hombres y mujeres).

El valor observado es 15 y 19. Para el calculo del valor esperado se obtiene de (15+19)/2=17. Resumiendose la informacion queda de la siguiente manera:

Solucion

frec<-c(15,19)
chisq.test(frec)

## 
##  Chi-squared test for given probabilities
## 
## data:  frec
## X-squared = 0.47059, df = 1, p-value = 0.4927

Calculo de los valores esperados

chisq.test(frec)$expected

## [1] 17 17

Calculo de los valores criticos

qchisq(0.95,1)

## [1] 3.841459

TEST DE INDEPENDENCIA - FISHERS TEST

Es una prueba de significacion estadistica, utilizadas en el analisis de tablas de contingencia cuando los tamalos de las muestras son pequenos (menores a 5)

matrix<-matrix(c(4,11,10,13,3,4,6,8),nrow=2)
matrix

##      [,1] [,2] [,3] [,4]
## [1,]    4   10    3    6
## [2,]   11   13    4    8

chisq.test(matrix)

## Warning in chisq.test(matrix): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  matrix
## X-squared = 1.2845, df = 3, p-value = 0.7328

Test de Fisher

fisher.test(matrix)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix
## p-value = 0.7229
## alternative hypothesis: two.sided

Como se observa los resultados de la prueba, resulta no significativa (p=0.7328). Asi mismo tenemos grados de libertad menores a 5. Cuando esto sucede es recomendable utilizar un test de independencia denominado Fishers test.

Como estamos trabajando una matriz de datos, es posible calcular las proporciones que corresponden los datos, tanto para las columnas como para las filas.

matrix

##      [,1] [,2] [,3] [,4]
## [1,]    4   10    3    6
## [2,]   11   13    4    8

prop.table(matrix) 

##            [,1]      [,2]       [,3]      [,4]
## [1,] 0.06779661 0.1694915 0.05084746 0.1016949
## [2,] 0.18644068 0.2203390 0.06779661 0.1355932

Proporcion de filas

prop.table(matrix, 1) # proporcion de las filas  

##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.1739130 0.4347826 0.1304348 0.2608696
## [2,] 0.3055556 0.3611111 0.1111111 0.2222222

Proporicion de las columnas

prop.table(matrix, 2) # proporcion de las columnas

##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.2666667 0.4347826 0.4285714 0.4285714
## [2,] 0.7333333 0.5652174 0.5714286 0.5714286

Ejercicios 1.

library(MASS) # Carga el paquete MASS
data(survey)
names(survey) # visualiza titulos de las columnas de los datos

##  [1] "Sex"    "Wr.Hnd" "NW.Hnd" "W.Hnd"  "Fold"   "Pulse"  "Clap"  
##  [8] "Exer"   "Smoke"  "Height" "M.I"    "Age"

cuadro<- table(survey$Smoke, survey$Sex) 
cuadro

##        
##         Female Male
##   Heavy      5    6
##   Never     99   89
##   Occas      9   10
##   Regul      5   12

Cuantos grados de libertad tengo?

Es posible que los datos sean independientes?

Forma grafica

mosaicplot(cuadro, color=TRUE, main="Plot de mosaico")

Ejemplo 2

Supongamos que se quiere estudiar la posible asociacion entre el hecho de que una gestante fume durante el embarazo y que el nino presente bajo peso al nacer. Por lo tanto, se trata de ver si la probabilidad de tener bajo peso es diferente en gestantes que fumen o en gestantes que no fumen durante la gestacion. Para responder a esta pregunta se realiza un estudio de seguimiento sobre una cohorte de 2000 gestantes, a las que se interroga sobre su habito durante la gestacion y se determina ademas el peso del recien nacido. Los resultados de este estudio se muestran a continuacion:

Cuadro. 1. Tabla de contingencia para estudiar la asociacion entre fumar durante la gestacion y el bajo peso del nino al nacer. Estudio de seguimiento de 2000 gestantes

Planteamiento de la hipotesis:

H0: No hay asociacion entre las variables (en el ejemplo, el bajo peso del nino y el hecho de fumar durante la gestacion son independientes o no estan asociados).

H1: Si hay asociacion entre las variables, es decir, el bajo peso y el fumar durante la gestacion estan asociados.

Si<-c(43,105)
No<-c(207,1645)
cuadro3<-data.frame(Si,No)
rownames(cuadro3)<-c("Fumadores","No fumadores")
cuadro3

##               Si   No
## Fumadores     43  207
## No fumadores 105 1645

chisq.test(cuadro3)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  cuadro3
## X-squared = 38.427, df = 1, p-value = 5.685e-10

Denote que como solo tenemos un grado de libertad, de forma automatica se realiza la correccion de Yates (util para cuando tablas de 2x2 (dos columnas por dos filas)). La correccion de Yates solo funciona cuando se tiene una tabla de contingencia.

fisher.test(cuadro3)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  cuadro3
## p-value = 1.687e-08
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.161142 4.830833
## sample estimates:
## odds ratio 
##    3.25158

mosaicplot(cuadro3, color=TRUE, main="Plot de mosaico")

Conclusion: S? hay asociacion entre las variables, es decir, el bajo peso y el fumar durante la gestacion estan asociados de manera significativa (chi=38.42, gl = 1, p-value <0.05)

ANDEVA (Analisis de varianza)

Una de las herramientas mas versatiles y utiles en bioestadiistica es el anaalisis de varianza, abreviado como ANDEVA.

ANDEVA de una via

Ejemplo 3

Numero de granos de polen en las flores de C. rupestris segun los tratamientos de adicion de nectar. Aumentado: adiciono nectar. Disminuido: Elimino todo el nectar. Control: flores sin manipulacion

aum <- c(42, 48, 55, 42, 55, 38, 42, 49, 44, 54)
dis <- c(26, 29, 27, 32, 30, 28, 29, 28, 32, 31)
ctl <- c(39, 34, 39, 40, 38, 36, 36, 40, 39, 37)
polen <- c(aum,dis,ctl)
polen

##  [1] 42 48 55 42 55 38 42 49 44 54 26 29 27 32 30 28 29 28 32 31 39 34 39
## [24] 40 38 36 36 40 39 37

trats <- gl(3,10,label=c("aumentado","disminuido","control"))
trats

##  [1] aumentado  aumentado  aumentado  aumentado  aumentado  aumentado 
##  [7] aumentado  aumentado  aumentado  aumentado  disminuido disminuido
## [13] disminuido disminuido disminuido disminuido disminuido disminuido
## [19] disminuido disminuido control    control    control    control   
## [25] control    control    control    control    control    control   
## Levels: aumentado disminuido control

length(trats)==length(polen)

## [1] TRUE

is.factor(trats)

## [1] TRUE

modelo.andeva <- aov(polen~trats)
summary(modelo.andeva)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## trats        2 1566.9   783.4   50.35 7.76e-10 ***
## Residuals   27  420.1    15.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 boxplot(polen~trats, col="gray")

Pruebas post hoc

TukeyHSD(modelo.andeva)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = polen ~ trats)
## 
## $trats
##                       diff        lwr        upr     p adj
## disminuido-aumentado -17.7 -22.073802 -13.326198 0.0000000
## control-aumentado     -9.1 -13.473802  -4.726198 0.0000576
## control-disminuido     8.6   4.226198  12.973802 0.0001228

plot(TukeyHSD(modelo.andeva))

Ejercicio 2.

Contiene los datos de un experimento en el que los participantes recibian una entre tres posibles dosis de una droga experimental, y tras la cual se media el grado de vigilancia en una tarea.

Vigilancia<-c(30,38,35,41,27,24,32,26,31,29,27,35,21,25,17,21,20,19)
Dosis<-gl(3,6,label=c("A","B","C"))
Dosis

##  [1] A A A A A A B B B B B B C C C C C C
## Levels: A B C

is.factor(Dosis)

## [1] TRUE

data.frame(Vigilancia,Dosis)

##    Vigilancia Dosis
## 1          30     A
## 2          38     A
## 3          35     A
## 4          41     A
## 5          27     A
## 6          24     A
## 7          32     B
## 8          26     B
## 9          31     B
## 10         29     B
## 11         27     B
## 12         35     B
## 13         21     C
## 14         25     C
## 15         17     C
## 16         21     C
## 17         20     C
## 18         19     C

anova1<-aov(Vigilancia~Dosis)
summary(anova1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Dosis        2    481   240.5   11.68 0.000876 ***
## Residuals   15    309    20.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Recuerde: El primer argumento es siempre la variable dependiente (Vigilancia), que es seguido por el simbolo (~) y despues la variable independiente.

Estadisticos descriptivos

tapply(Vigilancia,Dosis, mean)

##    A    B    C 
## 32.5 30.0 20.5

tapply(Vigilancia,Dosis, sd)

##        A        B        C 
## 6.595453 3.346640 2.664583

boxplot(Vigilancia~Dosis, col="gray")

Post-hoc

TukeyHSD(anova1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Vigilancia ~ Dosis)
## 
## $Dosis
##      diff        lwr       upr     p adj
## B-A  -2.5  -9.306495  4.306495 0.6158245
## C-A -12.0 -18.806495 -5.193505 0.0009912
## C-B  -9.5 -16.306495 -2.693505 0.0066219

plot(TukeyHSD(anova1))

ANOVA Factorial

Se trata de un disenpo de bloques aleatorizados (cada finca es un bloque). Introducimos los datos con las instrucciones

produc <- c(2.1, 2.2, 1.8, 2, 1.9, 2.2, 2.6, 2.7, 2.5, 2.8, 1.8, 1.9, 1.6, 2, 1.9, 2.1, 2, 2.2, 2.4, 2.1)
fert <- gl(4, 5)
finca <- factor(rep(1:5, 4))
xtabs(produc ~ finca + fert)

##      fert
## finca   1   2   3   4
##     1 2.1 2.2 1.8 2.1
##     2 2.2 2.6 1.9 2.0
##     3 1.8 2.7 1.6 2.2
##     4 2.0 2.5 2.0 2.4
##     5 1.9 2.8 1.9 2.1

tapply(produc, fert, summary)

## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.8     1.9     2.0     2.0     2.1     2.2 
## 
## $`2`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.20    2.50    2.60    2.56    2.70    2.80 
## 
## $`3`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.60    1.80    1.90    1.84    1.90    2.00 
## 
## $`4`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    2.10    2.10    2.16    2.20    2.40

tapply(produc, finca, summary)

## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.800   2.025   2.100   2.050   2.125   2.200 
## 
## $`2`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.900   1.975   2.100   2.175   2.300   2.600 
## 
## $`3`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.600   1.750   2.000   2.075   2.325   2.700 
## 
## $`4`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.000   2.000   2.200   2.225   2.425   2.500 
## 
## $`5`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.900   1.900   2.000   2.175   2.275   2.800

stripchart(produc ~ fert, method = "stack")

stripchart(produc ~ finca, method = "stack")

interaction.plot(fert, finca, produc, legend = T)

interaction.plot(finca, fert, produc, legend = T)

aovfinca <- lm(produc ~ finca + fert)
summary(aovfinca)

## 
## Call:
## lm(formula = produc ~ finca + fert)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2700 -0.1350  0.0250  0.1175  0.2050 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9100     0.1166  16.378 1.42e-09 ***
## finca2        0.1250     0.1304   0.959 0.356625    
## finca3        0.0250     0.1304   0.192 0.851151    
## finca4        0.1750     0.1304   1.342 0.204377    
## finca5        0.1250     0.1304   0.959 0.356625    
## fert2         0.5600     0.1166   4.802 0.000432 ***
## fert3        -0.1600     0.1166  -1.372 0.195171    
## fert4         0.1600     0.1166   1.372 0.195171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1844 on 12 degrees of freedom
## Multiple R-squared:  0.7884, Adjusted R-squared:  0.6649 
## F-statistic: 6.387 on 7 and 12 DF,  p-value: 0.00272

p.aov <- aov(produc ~ finca + fert)
efectos <- model.tables(p.aov)
efectos

## Tables of effects
## 
##  finca 
## finca
##      1      2      3      4      5 
## -0.090  0.035 -0.065  0.085  0.035 
## 
##  fert 
## fert
##     1     2     3     4 
## -0.14  0.42 -0.30  0.02

Ejercicio 3.

data(InsectSprays)
str(InsectSprays)

## 'data.frame':    72 obs. of  2 variables:
##  $ count: num  10 7 20 14 14 12 10 23 17 20 ...
##  $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

attach(InsectSprays)
tapply(count, spray, mean)

##         A         B         C         D         E         F 
## 14.500000 15.333333  2.083333  4.916667  3.500000 16.666667

tapply(count, spray, var)

##         A         B         C         D         E         F 
## 22.272727 18.242424  3.901515  6.265152  3.000000 38.606061

tapply(count, spray, length)

##  A  B  C  D  E  F 
## 12 12 12 12 12 12

boxplot(count ~ spray)

oneway.test(count ~ spray)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  count and spray
## F = 36.065, num df = 5.000, denom df = 30.043, p-value = 7.999e-12

aov.out = aov(count ~ spray, data=InsectSprays)
TukeyHSD(aov.out)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = count ~ spray, data = InsectSprays)
## 
## $spray
##            diff        lwr       upr     p adj
## B-A   0.8333333  -3.866075  5.532742 0.9951810
## C-A -12.4166667 -17.116075 -7.717258 0.0000000
## D-A  -9.5833333 -14.282742 -4.883925 0.0000014
## E-A -11.0000000 -15.699409 -6.300591 0.0000000
## F-A   2.1666667  -2.532742  6.866075 0.7542147
## C-B -13.2500000 -17.949409 -8.550591 0.0000000
## D-B -10.4166667 -15.116075 -5.717258 0.0000002
## E-B -11.8333333 -16.532742 -7.133925 0.0000000
## F-B   1.3333333  -3.366075  6.032742 0.9603075
## D-C   2.8333333  -1.866075  7.532742 0.4920707
## E-C   1.4166667  -3.282742  6.116075 0.9488669
## F-C  14.5833333   9.883925 19.282742 0.0000000
## E-D  -1.4166667  -6.116075  3.282742 0.9488669
## F-D  11.7500000   7.050591 16.449409 0.0000000
## F-E  13.1666667   8.467258 17.866075 0.0000000

plot(TukeyHSD(aov.out))

Comparacion de modelos

aov.out2<-aov(count ~1, data=InsectSprays)
summary(aov.out2)

##             Df Sum Sq Mean Sq F value Pr(>F)
## Residuals   71   3684   51.89

library(bbmle)

## Loading required package: stats4

AICctab(aov.out, aov.out2,  base= T, delta = T, sort = T, weights = T, nobs =56)

##          AICc  dAICc df weight
## aov.out  411.2   0.0 7  1     
## aov.out2 491.9  80.7 2  <0.001

library(visreg)
visreg(aov.out,"spray", partial = F)

library(car)
par(mfrow=c(2,2))
lm.r<-lm(count ~ spray, data=InsectSprays)
plot(lm.r)

Supuestos del Modelo ANOVA

residuals(lm.r)->residuos
shapiro.test(residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96006, p-value = 0.02226

plot(predict(lm.r), residuals(lm.r)) #grafica los residuos
abline(h=0,lty=2,col="red")

ncvTest(lm.r)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 20.00064    Df = 1     p = 7.741633e-06

library(car)
levene.test(count ~ spray, data=InsectSprays)

## Warning: 'levene.test' is deprecated.
## Use 'leveneTest' instead.
## See help("Deprecated") and help("car-deprecated").

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  5  3.8214 0.004223 **
##       66                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No parametria

kruskal.test(count ~ spray, data=InsectSprays)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  count by spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10

crec <- c(4.39, 4.40, 2.96, 4.58, 4.83, 3.64, 4.98, 2.90, 4.47, 3.44, 6.96, 5.49, 6.04,
          6.19, 6.73, 4.09, 5.53, 6.59, 7.04, 5.63, 3.92, 5.03, 4.20, 5.25, 4.22, 5.28, 4.42, 4.43,
          4.08, 4.81, 4.51, 5.85, 5.30, 4.07, 6.21, 6.34, 6.54, 4.49, 6.39, 3.75)
luz <- gl(2,20,labels=c("Sol", "Sombra"))
altura <- gl(2, 10, length=40, labels=c("Sub", "Paramo"))
cal <- data.frame(crec,luz,altura)
cal

##    crec    luz altura
## 1  4.39    Sol    Sub
## 2  4.40    Sol    Sub
## 3  2.96    Sol    Sub
## 4  4.58    Sol    Sub
## 5  4.83    Sol    Sub
## 6  3.64    Sol    Sub
## 7  4.98    Sol    Sub
## 8  2.90    Sol    Sub
## 9  4.47    Sol    Sub
## 10 3.44    Sol    Sub
## 11 6.96    Sol Paramo
## 12 5.49    Sol Paramo
## 13 6.04    Sol Paramo
## 14 6.19    Sol Paramo
## 15 6.73    Sol Paramo
## 16 4.09    Sol Paramo
## 17 5.53    Sol Paramo
## 18 6.59    Sol Paramo
## 19 7.04    Sol Paramo
## 20 5.63    Sol Paramo
## 21 3.92 Sombra    Sub
## 22 5.03 Sombra    Sub
## 23 4.20 Sombra    Sub
## 24 5.25 Sombra    Sub
## 25 4.22 Sombra    Sub
## 26 5.28 Sombra    Sub
## 27 4.42 Sombra    Sub
## 28 4.43 Sombra    Sub
## 29 4.08 Sombra    Sub
## 30 4.81 Sombra    Sub
## 31 4.51 Sombra Paramo
## 32 5.85 Sombra Paramo
## 33 5.30 Sombra Paramo
## 34 4.07 Sombra Paramo
## 35 6.21 Sombra Paramo
## 36 6.34 Sombra Paramo
## 37 6.54 Sombra Paramo
## 38 4.49 Sombra Paramo
## 39 6.39 Sombra Paramo
## 40 3.75 Sombra Paramo

replications(crec~luz*altura, data=cal)

##        luz     altura luz:altura 
##         20         20         10

is.factor(altura)

## [1] TRUE

m1 <- aov(crec~luz*altura, data=cal)
summary(m1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## luz          1  0.080   0.080   0.117   0.7346    
## altura       1 18.920  18.920  27.564 6.98e-06 ***
## luz:altura   1  3.534   3.534   5.149   0.0294 *  
## Residuals   36 24.711   0.686                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model.tables(m1, "means")

## Tables of means
## Grand mean
##         
## 4.99925 
## 
##  luz 
## luz
##    Sol Sombra 
##  5.044  4.954 
## 
##  altura 
## altura
##    Sub Paramo 
##  4.311  5.687 
## 
##  luz:altura 
##         altura
## luz      Sub   Paramo
##   Sol    4.059 6.029 
##   Sombra 4.564 5.345

model.tables(m1, "effects")

## Tables of effects
## 
##  luz 
## luz
##      Sol   Sombra 
##  0.04475 -0.04475 
## 
##  altura 
## altura
##     Sub  Paramo 
## -0.6877  0.6877 
## 
##  luz:altura 
##         altura
## luz      Sub      Paramo  
##   Sol    -0.29725  0.29725
##   Sombra  0.29725 -0.29725

interaction.plot(cal$altura,cal$luz,cal$crec)

with(cal,interaction.plot(altura,luz,crec))

boxplot(crec~luz*altura, data=cal, col="forestgreen")

GLM

library(AICcmodavg)

## 
## Attaching package: 'AICcmodavg'
## 
## The following object is masked from 'package:bbmle':
## 
##     AICc

library(ggplot2)

library(AICcmodavg)
library(ggplot2)
data(lizards)
head(lizards)

##   Insolation Diameter Height      Time  Species Counts
## 1      sunny      <=2     <5   morning  grahami     20
## 2      sunny      <=2     <5   morning opalinus      2
## 3      sunny      <=2     <5    midday  grahami      8
## 4      sunny      <=2     <5    midday opalinus      1
## 5      sunny      <=2     <5 afternoon  grahami      4
## 6      sunny      <=2     <5 afternoon opalinus      4

GLM.1<- glm(Counts ~ Insolation + Diameter + Height, family = "poisson" , data = lizards)
summary(GLM.1)

## 
## Call:
## glm(formula = Counts ~ Insolation + Diameter + Height, family = "poisson", 
##     data = lizards)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.3007  -2.5045  -1.0987   0.8778   7.4309  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.41280    0.06532  52.251  < 2e-16 ***
## Insolationsunny -1.48684    0.10858 -13.694  < 2e-16 ***
## Diameter>2      -0.45447    0.08640  -5.260 1.44e-07 ***
## Height>=5       -0.60659    0.08812  -6.884 5.82e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 737.56  on 47  degrees of freedom
## Residual deviance: 416.89  on 44  degrees of freedom
## AIC: 588.14
## 
## Number of Fisher Scoring iterations: 5

plot(GLM.1)

exp(-1.4868)#Insolationsunny

## [1] 0.226095

exp(-0.45447) #Diameter>2

## [1] 0.6347843

exp(-0.60659 ) #Height>=5

## [1] 0.5452069

1-pchisq(737.56- 416.89, df=(47-44)) #Logistic regression model provides an adequate fit for the data

## [1] 0

De cada lagartija que este en un lugar sombreado, disminuye en 0.22 las lagartijas que estan en un lugar soleado.

GLM.2<- glm(Counts ~ Insolation * Diameter* Height, , family = "poisson" , data = lizards)
summary(GLM.2)

## 
## Call:
## glm(formula = Counts ~ Insolation * Diameter * Height, family = "poisson", 
##     data = lizards)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.482  -2.226  -1.218   1.316   6.845  
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)                           3.29584    0.07857  41.949  < 2e-16
## Insolationsunny                      -1.42403    0.17836  -7.984 1.42e-15
## Diameter>2                           -0.14595    0.11539  -1.265 0.205931
## Height>=5                            -0.37807    0.12321  -3.068 0.002152
## Insolationsunny:Diameter>2           -0.33955    0.28394  -1.196 0.231753
## Insolationsunny:Height>=5             0.21101    0.26669   0.791 0.428820
## Diameter>2:Height>=5                 -0.71343    0.20881  -3.417 0.000634
## Insolationsunny:Diameter>2:Height>=5 -0.21813    0.51596  -0.423 0.672466
##                                         
## (Intercept)                          ***
## Insolationsunny                      ***
## Diameter>2                              
## Height>=5                            ** 
## Insolationsunny:Diameter>2              
## Insolationsunny:Height>=5               
## Diameter>2:Height>=5                 ***
## Insolationsunny:Diameter>2:Height>=5    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 737.56  on 47  degrees of freedom
## Residual deviance: 396.07  on 40  degrees of freedom
## AIC: 575.32
## 
## Number of Fisher Scoring iterations: 6

Regresion

grasas <- read.table("http://www.uam.es/joser.berrendero/datos/EdadPesoGrasas.txt",  header = TRUE)
names(grasas)

## [1] "peso"   "edad"   "grasas"

pairs(grasas)

cor(grasas)

##             peso      edad    grasas
## peso   1.0000000 0.2400133 0.2652935
## edad   0.2400133 1.0000000 0.8373534
## grasas 0.2652935 0.8373534 1.0000000

regresion <- lm(grasas ~ edad, data = grasas)
regresion2 <- lm(grasas ~ edad+peso, data = grasas)
summary(regresion)

## 
## Call:
## lm(formula = grasas ~ edad, data = grasas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.478 -26.816  -3.854  28.315  90.881 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 102.5751    29.6376   3.461  0.00212 ** 
## edad          5.3207     0.7243   7.346 1.79e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.46 on 23 degrees of freedom
## Multiple R-squared:  0.7012, Adjusted R-squared:  0.6882 
## F-statistic: 53.96 on 1 and 23 DF,  p-value: 1.794e-07

Coeficiente de determinacion. Mide variabilidad explicada de Y sobre X

anova(regresion)

## Analysis of Variance Table
## 
## Response: grasas
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## edad       1 101933  101933  53.964 1.794e-07 ***
## Residuals 23  43444    1889                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mide la bondad del ajuste de la recta a los datos.

plot(grasas$edad, grasas$grasas, xlab = "Edad", ylab = "Grasas")
abline(regresion)

par(mfrow=c(2,2))
plot(lm(grasas ~ edad, data = grasas))

Residual vs Fitted: Deberia estar distribuidos aleatoriamente alrededor de la linea horizontal que representa un error residual de cero

Normal Q-Q: debera sugerir que los errores residuales se distribuyen normalmente

Scale-Location Muestra la raiz cuadrada de los residuos estandarizados, como una funcion de los valores ajustados. debera haber una tendencia clara en ese trama

Residual vs Leverage Las distancias mas grandes que 1 son sospechosos y sugieren la presencia de un valor at??pico posible. y su eliminacion podria tener efectos sobre la regresion

Queremos utilizar la recta de minimos cuadrados para predecir la cantidad de grasas para individuos de edades 31,31,32,33,34,35

nuevas.edades <- data.frame(edad = seq(30, 35))
nuevas.edades

##   edad
## 1   30
## 2   31
## 3   32
## 4   33
## 5   34
## 6   35

predict(regresion, nuevas.edades)

##        1        2        3        4        5        6 
## 262.1954 267.5161 272.8368 278.1575 283.4781 288.7988

predict(regresion,data.frame(nuevas.edades), level = 0.95,interval = "confidence") 

##        fit      lwr      upr
## 1 262.1954 239.6112 284.7797
## 2 267.5161 245.8056 289.2266
## 3 272.8368 251.9291 293.7445
## 4 278.1575 257.9731 298.3419
## 5 283.4781 263.9288 303.0275
## 6 288.7988 269.7874 307.8102

Por ejemplo, para un individuo de 30 anos, predecimos una cantidad de grasas de 262.1954

confint(regresion,level=0.95)

##                 2.5 %     97.5 %
## (Intercept) 41.265155 163.885130
## edad         3.822367   6.818986

Otros analisis

Outliers-QQ-Plots

library(car)
outlierTest(regresion) #Outliers???Bonferonnitest 

## 
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferonni p
## 8 2.334884           0.029076       0.7269

influenceIndexPlot(regresion,id.n=3) #High leverage (hat) points (graph)

influencePlot(regresion,id.n=3)

##       StudRes        Hat     CookD
## 2  -0.4634322 0.14153039 0.1353869
## 5   1.1189868 0.12878822 0.3025650
## 6   1.6675344 0.09537193 0.3688438
## 8   2.3348835 0.04270352 0.3191819
## 15 -0.6654287 0.16108247 0.2087260
## 16  1.9069349 0.06190011 0.3280782

qqPlot(regresion)#testting fornormality 

ncvTest(regresion) #Prueba de heterocedasticidad, Null is constant variance 

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.04653681    Df = 1     p = 0.8292029

(Se prueba si,varianza de las perturbaciones no es constante a lo largo de las observaciones.)

residuos <- rstandard(regresion)
valores.ajustados <- fitted(regresion)
plot(valores.ajustados, residuos)

No se observa ningun patron especial, por lo que tanto

la homocedasticidad como la linealidad resultan hipotesis razonables.

qqnorm(residuos)
qqline(residuos)

plot(regresion)