# R-code for Lecture 8, vt19. # 6/5-19. AL #################### Fruit data again: fruit <- data.frame(type = c(rep(0, 300), rep(1, 70), rep(2, 40)), green = c(rep(1, 100), rep(0, 200), rep(1, 40), rep(0, 30), rep(1, 20), rep(0, 20))) fruit\$type <- factor(fruit\$type, levels = c(0, 1, 2), labels = c("banana", "apple", "melon")) fruit\$green <- factor(fruit\$green, levels = c(0, 1), labels = c("yellow", "green")) n <- nrow(fruit) # Fit the model logit(p_green) = b0 + b1*apple + b2*melon model <- glm(green ~ type, data = fruit, family = "binomial") summary(model) #Likelihood ratio test: (Ddiff <- model\$null.deviance - model\$deviance) (dfdiff <- model\$df.null - model\$df.residual) # .. or using the deviance table anova(update(model, . ~ 1), model) # chi2-quantile for comparison qchisq(1 - 0.05, dfdiff) # P-value for the test: pchisq(Ddiff, dfdiff, lower.tail = FALSE) # Comparing using AIC/BIC # Are all three fruit tytes different from each other, or is just one of them # different from the other two? Take the null model and add indicators for # each separate fruit type: model.0 <- update(model, . ~ 1) model.banana <- update(model.0, . ~ I(type != "banana")) model.apple <- update(model.0, . ~ I(type != "apple")) model.melon <- update(model.0, . ~ I(type != "melon")) model.full <- model cbind(AIC(model.0, model.banana, model.apple, model.melon, model.full), BIC(model.0, model.banana, model.apple, model.melon, model.full)) # Pseudo R2 # Likelihood = exp(-deviance/2) L0 <- exp(-model.0\$deviance/2) R2max <- 1 - L0^(2/n) L <- exp(-c(model.0\$deviance, model.banana\$deviance, model.apple\$deviance, model.melon\$deviance, model\$deviance)/2) R2CS <- 1 - (L0/L)^(2/n) R2N <- R2CS/R2max R2 <- round(100*cbind(R2CS, R2N), 2) row.names(R2) <- c("model.0", "model.banana", "model.apple", "model.melon", "model.full") R2 100*R2max