#### GENERAL LINEAR MODEL > rabbits1 <- readXL("D:/Nina/statistica/stat 2019-2020/занятие 8/linear.models.xlsx", rownames=FALSE, header=TRUE, na="", sheet="logit", stringsAsFactors=TRUE) > summary (rabbits1) femaleID litterID weight fat parasites density disturbance color survival habitat n1003 : 1 n37 : 4 Min. :115.0 Min. : 4.500 Min. :2.000 Min. : 12.12 Min. : 20.12 black:17 Min. :0.0000 a:28 n1048 : 1 n135 : 3 1st Qu.:153.0 1st Qu.: 7.500 1st Qu.:4.000 1st Qu.: 30.68 1st Qu.: 87.46 grey :18 1st Qu.:0.0000 b:27 n1055 : 1 n15 : 3 Median :188.0 Median : 8.659 Median :5.000 Median : 36.45 Median :206.48 white:20 Median :1.0000 n1086 : 1 n2 : 3 Mean :188.6 Mean : 8.977 Mean :5.127 Mean : 41.43 Mean :191.57 Mean :0.6364 n1164 : 1 n28 : 3 3rd Qu.:214.5 3rd Qu.:10.360 3rd Qu.:6.000 3rd Qu.: 49.54 3rd Qu.:303.52 3rd Qu.:1.0000 n1205 : 1 n33 : 3 Max. :319.0 Max. :13.625 Max. :9.000 Max. :127.60 Max. :384.01 Max. :1.0000 (Other):49 (Other):36 > rabbit1 <- lm(weight ~ fat + color + habitat + color:habitat, data = rabbits1) > summary (rabbit1) Call: lm(formula = weight ~ fat + color + habitat + color:habitat, data = rabbits1) Residuals: Min 1Q Median 3Q Max -84.199 -22.016 -4.883 25.369 112.873 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 164.787 30.597 5.386 0.00000214 *** fat 4.174 2.881 1.449 0.15388 color[T.grey] -28.900 18.506 -1.562 0.12494 color[T.white] -37.497 18.875 -1.987 0.05270 . habitat[T.b] -34.653 18.958 -1.828 0.07379 . color[T.grey]:habitat[T.b] 74.183 26.369 2.813 0.00709 ** color[T.white]:habitat[T.b] 78.921 28.400 2.779 0.00776 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 38.63 on 48 degrees of freedom Multiple R-squared: 0.3229, Adjusted R-squared: 0.2383 F-statistic: 3.815 on 6 and 48 DF, p-value: 0.003451 > anova(rabbit1) Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) fat 1 14926 14925.8 10.0014 0.002711 ** color 2 926 462.9 0.3102 0.734785 habitat 1 3023 3023.1 2.0257 0.161124 color:habitat 2 15288 7644.0 5.1221 0.009632 ** Residuals 48 71634 1492.4 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > plot (rabbit1) ############ BEST MODEL SELECTION > library(MuMIn, pos=22) > options(na.action = na.fail) > res1 <- dredge(rabbit1, rank="AICc", REML=F) > est1 <- model.avg(res1, subset=delta<7, revised.var=T) > summary(est1) Call: model.avg(object = res1, subset = delta < 7, revised.var = T) Component model call: lm(formula = weight ~ <8 unique rhs>, data = rabbits1) Component models: df logLik AICc delta weight 134 7 -276.45 569.28 0.00 0.31 1234 8 -275.27 569.67 0.39 0.25 2 3 -281.81 570.10 0.82 0.20 23 4 -280.84 570.47 1.19 0.17 12 5 -281.53 574.29 5.01 0.02 3 3 -284.08 574.63 5.35 0.02 123 6 -280.59 574.93 5.65 0.02 (Null) 2 -286.00 576.22 6.94 0.01 Term codes: color fat habitat color:habitat 1 2 3 4 Model-averaged coefficients: (full average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 158.284 43.528 43.837 3.611 0.000305 *** color[T.grey] -17.121 21.525 21.761 0.787 0.431411 color[T.white] -23.623 25.652 25.849 0.914 0.360774 habitat[T.b] -17.036 27.318 27.530 0.619 0.536051 color[T.grey]:habitat[T.b] 42.699 43.004 43.237 0.988 0.323370 color[T.white]:habitat[T.b] 49.359 48.980 49.194 1.003 0.315694 fat 4.116 3.915 3.947 1.043 0.297012 (conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 158.284 43.528 43.837 3.611 0.000305 *** color[T.grey] -28.603 21.109 21.510 1.330 0.183596 color[T.white] -39.466 21.774 22.159 1.781 0.074904 . habitat[T.b] -22.334 29.326 29.586 0.755 0.450322 color[T.grey]:habitat[T.b] 76.858 26.522 27.194 2.826 0.004709 ** color[T.white]:habitat[T.b] 88.846 28.461 29.117 3.051 0.002278 ** fat 6.194 3.193 3.251 1.905 0.056747 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Relative variable importance: habitat fat color color:habitat Importance: 0.76 0.66 0.60 0.56 N containing models: 5 5 4 2 ############## LIKELIHOOD RATIO TESTS ## FAT EFFECT > library(foreign, pos=23) > library(survival, pos=23) > library(MASS, pos=23) > library(epiDisplay, pos=23) > rabbit2 <- lm(weight ~ color + habitat + color:habitat, data = rabbits1) > lrtest (rabbit1, rabbit2) Likelihood ratio test Model 1: weight ~ fat + color + habitat + color:habitat Model 2: weight ~ color + habitat + color:habitat #Df LogLik Df Chisq Pr(>Chisq) 1 8 -275.27 2 7 -276.45 -1 2.3541 0.125 > rabbit3 <- lm(weight ~ fat + color + habitat, data = rabbits1) > lrtest (rabbit1, rabbit3) Likelihood ratio test Model 1: weight ~ fat + color + habitat + color:habitat Model 2: weight ~ fat + color + habitat #Df LogLik Df Chisq Pr(>Chisq) 1 8 -275.27 2 6 -280.59 -2 10.639 0.004894 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ####HIGHTLY SIGNIFICANT > library(emmeans, pos=27) > lsmeans(rabbit1, pairwise~color*habitat, adjust="tukey") $lsmeans color habitat lsmean SE df lower.CL upper.CL black a 202 13.0 48 176 228 grey a 173 13.0 48 147 199 white a 165 13.0 48 139 191 black b 168 13.7 48 140 195 grey b 213 12.9 48 187 239 white b 209 13.3 48 182 236 Confidence level used: 0.95 $contrasts contrast estimate SE df t.ratio p.value black,a - grey,a 28.90 18.5 48 1.562 0.6270 black,a - white,a 37.50 18.9 48 1.987 0.3649 black,a - black,b 34.65 19.0 48 1.828 0.4582 black,a - grey,b -10.63 18.4 48 -0.578 0.9920 black,a - white,b -6.77 18.1 48 -0.375 0.9990 grey,a - white,a 8.60 18.0 48 0.477 0.9967 grey,a - black,b 5.75 18.8 48 0.306 0.9996 grey,a - grey,b -39.53 18.2 48 -2.169 0.2710 grey,a - white,b -35.67 19.0 48 -1.880 0.4263 white,a - black,b -2.84 18.7 48 -0.152 1.0000 white,a - grey,b -48.13 18.2 48 -2.650 0.1049 white,a - white,b -44.27 19.9 48 -2.227 0.2445 black,b - grey,b -45.28 18.8 48 -2.412 0.1725 black,b - white,b -41.42 19.3 48 -2.146 0.2818 grey,b - white,b 3.86 18.7 48 0.206 0.9999 P value adjustment: tukey method for comparing a family of 6 estimates > with(rabbits1, plotMeans(weight, color, habitat, error.bars="se", connect=TRUE, legend.pos="farright")) ############################### ##########SCRIPT############ #### GENERAL LINEAR MODEL rabbits1 <- readXL("D:/Nina/statistica/stat 2019-2020/занятие 8/linear.models.xlsx", rownames=FALSE, header=TRUE, na="", sheet="logit", stringsAsFactors=TRUE) summary (rabbits1) rabbit1 <- lm(weight ~ fat + color + habitat + color:habitat, data = rabbits1) summary (rabbit1) anova(rabbit1) plot (rabbit1) ############ BEST MODEL SELECTION library(MuMIn, pos=22) options(na.action = na.fail) res1 <- dredge(rabbit1, rank="AICc", REML=F) est1 <- model.avg(res1, subset=delta<7, revised.var=T) summary(est1) ############## LIKELIHOOD RATIO TESTS ## FAT EFFECT library(foreign, pos=23) library(survival, pos=23) library(MASS, pos=23) library(epiDisplay, pos=23) rabbit2 <- lm(weight ~ color + habitat + color:habitat, data = rabbits1) lrtest (rabbit1, rabbit2) ## INTERACTION EFFECT rabbit3 <- lm(weight ~ fat + color + habitat, data = rabbits1) lrtest (rabbit1, rabbit3) ####HIGHTLY SIGNIFICANT #### TUKEY POST HOC COMPARISON library(emmeans, pos=27) lsmeans(rabbit1, pairwise~color*habitat, adjust="tukey") with(rabbits1, plotMeans(weight, color, habitat, error.bars="se", connect=TRUE, legend.pos="farright"))