packages

library(dplyr) 
library(ggpubr)
library(rcompanion)
library(car)

IPOMOPSIS.TXT

Este archivo se encuentra en un folder con datos de Crawley (2003) que adjunto aquí (Crawley). No hemos visto en clase este procedimiento para importar bases de datos. Para que este método funcione el folder adjunto debe de estar en el directorio de trabajo de RStudio (veremos esto en la siguiente sesión).

Leyendo los datos

read.table("Crawley/ipomopsis.txt", header=TRUE) -> ipo
names(ipo)
[1] "Root"    "Fruit"   "Grazing"
head(ipo)
tail(ipo)

PLOTS

plot(ipo$Fruit ~ ipo$Root, pch=20)

ggplot(aes(x = Grazing, y = Fruit, col= Grazing), data= ipo) +
geom_boxplot(outlier.colour = 2, orientation = "x") +
geom_point(position = position_dodge2(width = 0.2))

###
p <- ggplot(aes(x = Root, y = Fruit, color = Grazing), data= ipo)
p + geom_point() +
geom_smooth(method = "lm")
`geom_smooth()` using formula 'y ~ x'

p + geom_boxplot(outlier.colour = 2) +
geom_point() +
geom_smooth(method = "lm")
`geom_smooth()` using formula 'y ~ x'

Interacción entre “Root” y “Grazing” ?

ancova <- lm(ipo$Fruit ~ ipo$Grazing *
ipo$Root)
summary(ancova)

Call:
lm(formula = ipo$Fruit ~ ipo$Grazing * ipo$Root)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.3177  -2.8320   0.1247   3.8511  17.1313 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  -125.173     12.811  -9.771 1.15e-11 ***
ipo$GrazingUngrazed            30.806     16.842   1.829   0.0757 .  
ipo$Root                       23.240      1.531  15.182  < 2e-16 ***
ipo$GrazingUngrazed:ipo$Root    0.756      2.354   0.321   0.7500    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.831 on 36 degrees of freedom
Multiple R-squared:  0.9293,    Adjusted R-squared:  0.9234 
F-statistic: 157.6 on 3 and 36 DF,  p-value: < 2.2e-16

ANCOVA “species.txt”

read.table("Crawley/species.txt", header= TRUE) -> sp
names(sp)
[1] "pH"      "Biomass" "Species"
head(sp)
tail(sp)

Explorando datos

ggdensity(sp$Species)

shapiro.test(sp$Species)

    Shapiro-Wilk normality test

data:  sp$Species
W = 0.97805, p-value = 0.1317
p <- ggplot(aes(x = Biomass, y = Species, color = pH), data= sp)
p + geom_point() +
geom_smooth(method = "lm")

Modelo

ancova <- lm(sp$Species ~ sp$Biomass *
sp$pH)
summary(ancova)

Call:
lm(formula = sp$Species ~ sp$Biomass * sp$pH)

Residuals:
   Min     1Q Median     3Q    Max 
-9.290 -2.554 -0.124  2.208 15.677 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          40.60407    1.36701  29.703  < 2e-16 ***
sp$Biomass           -2.80045    0.23856 -11.739  < 2e-16 ***
sp$pHlow            -22.75667    1.83564 -12.397  < 2e-16 ***
sp$pHmid            -11.57307    1.86926  -6.191  2.1e-08 ***
sp$Biomass:sp$pHlow  -0.02733    0.51248  -0.053    0.958    
sp$Biomass:sp$pHmid   0.23535    0.38579   0.610    0.543    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.818 on 84 degrees of freedom
Multiple R-squared:  0.8531,    Adjusted R-squared:  0.8444 
F-statistic: 97.58 on 5 and 84 DF,  p-value: < 2.2e-16

PEERO ESTOS SON CONTEOS!!!!!

Mejor basarse en la Devianza no la varianza (“Poisson errors”)

p + geom_point() +
geom_smooth(method ='glm',method.args = list(family = 'poisson'))
`geom_smooth()` using formula 'y ~ x'

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQoNCiMjIyMgcGFja2FnZXMNCmBgYHtyLCBtZXNzYWdlPUZBTFNFfQ0KbGlicmFyeShkcGx5cikgDQpsaWJyYXJ5KGdncHVicikNCmxpYnJhcnkocmNvbXBhbmlvbikNCmxpYnJhcnkoY2FyKQ0KYGBgICANCg0KIyMjIyBJUE9NT1BTSVMuVFhUDQpFc3RlIGFyY2hpdm8gc2UgZW5jdWVudHJhIGVuIHVuIGZvbGRlciBjb24gZGF0b3MgZGUgQ3Jhd2xleSAoMjAwMykgcXVlIGFkanVudG8gYXF1w60gKFtDcmF3bGV5XShodHRwczovL2RyaXZlLmdvb2dsZS5jb20vZHJpdmUvZm9sZGVycy8xalJSWGRsaVJfRExzbTliV3p6SmNRUVNPeW1NR19Lckg/dXNwPXNoYXJpbmcpKS4gTm8gaGVtb3MgdmlzdG8gZW4gY2xhc2UgZXN0ZSBwcm9jZWRpbWllbnRvIHBhcmEgaW1wb3J0YXIgYmFzZXMgZGUgZGF0b3MuIFBhcmEgcXVlIGVzdGUgbcOpdG9kbyBmdW5jaW9uZSBlbCBmb2xkZXIgYWRqdW50byBkZWJlIGRlIGVzdGFyIGVuIGVsIGRpcmVjdG9yaW8gZGUgdHJhYmFqbyBkZSBSU3R1ZGlvICh2ZXJlbW9zIGVzdG8gZW4gbGEgc2lndWllbnRlIHNlc2nDs24pLg0KDQojIyMjIExleWVuZG8gbG9zIGRhdG9zDQpgYGB7cn0NCnJlYWQudGFibGUoIkNyYXdsZXkvaXBvbW9wc2lzLnR4dCIsIGhlYWRlcj1UUlVFKSAtPiBpcG8NCm5hbWVzKGlwbykNCmhlYWQoaXBvKQ0KdGFpbChpcG8pDQpgYGAgIA0KDQojIyMjIFBMT1RTDQpgYGB7ciwgbWVzc2FnZT1GQUxTRX0NCnBsb3QoaXBvJEZydWl0IH4gaXBvJFJvb3QsIHBjaD0yMCkNCmdncGxvdChhZXMoeCA9IEdyYXppbmcsIHkgPSBGcnVpdCwgY29sPSBHcmF6aW5nKSwgZGF0YT0gaXBvKSArDQpnZW9tX2JveHBsb3Qob3V0bGllci5jb2xvdXIgPSAyLCBvcmllbnRhdGlvbiA9ICJ4IikgKw0KZ2VvbV9wb2ludChwb3NpdGlvbiA9IHBvc2l0aW9uX2RvZGdlMih3aWR0aCA9IDAuMikpDQojIyMNCnAgPC0gZ2dwbG90KGFlcyh4ID0gUm9vdCwgeSA9IEZydWl0LCBjb2xvciA9IEdyYXppbmcpLCBkYXRhPSBpcG8pDQpwICsgZ2VvbV9wb2ludCgpICsNCmdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpDQpwICsgZ2VvbV9ib3hwbG90KG91dGxpZXIuY29sb3VyID0gMikgKw0KZ2VvbV9wb2ludCgpICsNCmdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpDQpgYGAgIA0KDQojIyMjIEludGVyYWNjacOzbiBlbnRyZSAiUm9vdCIgeSAiR3JhemluZyIgPw0KYGBge3J9DQphbmNvdmEgPC0gbG0oaXBvJEZydWl0IH4gaXBvJEdyYXppbmcgKg0KaXBvJFJvb3QpDQpzdW1tYXJ5KGFuY292YSkNCmBgYCAgDQoNCg0KIyMjIEFOQ09WQSAic3BlY2llcy50eHQiDQpgYGB7cn0NCnJlYWQudGFibGUoIkNyYXdsZXkvc3BlY2llcy50eHQiLCBoZWFkZXI9IFRSVUUpIC0+IHNwDQpuYW1lcyhzcCkNCmhlYWQoc3ApDQp0YWlsKHNwKQ0KYGBgICANCg0KIyMjIyBFeHBsb3JhbmRvIGRhdG9zDQpgYGB7cn0NCmdnZGVuc2l0eShzcCRTcGVjaWVzKQ0Kc2hhcGlyby50ZXN0KHNwJFNwZWNpZXMpDQpwIDwtIGdncGxvdChhZXMoeCA9IEJpb21hc3MsIHkgPSBTcGVjaWVzLCBjb2xvciA9IHBIKSwgZGF0YT0gc3ApDQpwICsgZ2VvbV9wb2ludCgpICsNCmdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpDQpgYGAgIA0KDQojIyMjIE1vZGVsbyANCmBgYHtyfQ0KYW5jb3ZhIDwtIGxtKHNwJFNwZWNpZXMgfiBzcCRCaW9tYXNzICoNCnNwJHBIKQ0Kc3VtbWFyeShhbmNvdmEpDQpgYGAgIA0KDQojIyMgUEVFUk8gRVNUT1MgU09OIENPTlRFT1MhISEhIQ0KIyMjIE1lam9yIGJhc2Fyc2UgZW4gbGEgRGV2aWFuemEgbm8gbGEgdmFyaWFuemEgKCJQb2lzc29uIGVycm9ycyIpDQoNCmBgYHtyfQ0KcCArIGdlb21fcG9pbnQoKSArDQpnZW9tX3Ntb290aChtZXRob2QgPSdnbG0nLG1ldGhvZC5hcmdzID0gbGlzdChmYW1pbHkgPSAncG9pc3NvbicpKQ0KYGBgICANCg0K