library(MASS) # Generierung der Daten mrdata set.seed(1234) x1 <- round(rnorm(20, 50, 5), 2) x2 <- round(rnorm(20, 200, 10), 2) x3 <- round(rnorm(20, 50, 5), 2) x4 <- round(rnorm(20, 100, 5), 2) X <- cbind(x1,x2,x3,x4) y <- numeric(20) for(i in 1:20){ y[i] <- 1000 - 20 * x1[i] + 5 * x2[i] + 15 * x3[i] + 1e-10 * x4[i] + rnorm(1, sd = 1) } y[4] <- min(y)+10 # vorher: 1986.426, nachher: 1410.371 y[8] <- max(y)-300 # vorher: 1659.097, nachher: 1686.426 data <- as.data.frame(cbind(X, y)) beta <- c(1000, -20, 5, 15, 1e-10) # (a) Explorative Untersuchung plot(data) color <- c(rep("black", 3), "red", rep("black", 3), "orange", rep("black",2), "green", rep("black", 5), "blue", rep("black", 2), "cyan") plot(data$x1,data$y, col=color, pch = 19, xlab=expression( x[1]), ylab =expression(y)) plot(data$x2,data$y, col=color, pch = 19, xlab=expression( x[2]), ylab =expression(y)) plot(data$x3,data$y, col=color, pch = 19, xlab=expression( x[3]), ylab =expression(y)) plot(data$x4,data$y, col=color, pch = 19, xlab=expression( x[4]), ylab =expression(y)) # (b) modelKQ <- lm(y ~ x1+x2+x3+x4, data) summary(modelKQ) round(abs(modelKQ$coefficients - beta)/beta, digits = 2) # Bemerkung: relativer Fehler für beta4 unsinnig, da fast bei 0 ! # (c) plot(1:20, modelKQ$residuals, pch = 16, col = color, xlab = " ", ylab = "Residuen", main = "Residuen-Plot fuer KQ-Schaetzer") abline(h = 0, lty = 2) for(i in 1:20){ lines(c(i,i), c(0,modelKQ$residuals[i])) } plot(1:20, abs(modelKQ$residuals), pch = 16, col = color, xlab = " ", ylab = "absolute Residuen", main = "Residuen-Plot fuer KQ-Schaetzer") abline(h = 0, lty = 2) for(i in 1:20){ lines(c(i,i), c(0,abs(modelKQ$residuals[i]))) } # (d) modelR <- rlm(y~x1+x2+x3+x4, data, psi = psi.hampel) # Standardeinstellungen der Tune-Parameter summary(modelR) round(abs(modelR$coefficients - beta)/beta, digits = 2) # Bemerkung: relativer Fehler für beta4 unsinnig, da fast bei 0 ! # (e) plot(1:20, modelR$residuals, pch = 16, col = color, xlab = " ", ylab = "Residuen", main = "Residuen-Plot fuer M-Schaetzer") abline(h = 0, lty = 2) for(i in 1:20){ lines(c(i,i), c(0,modelR$residuals[i])) } plot(1:20, abs(modelR$residuals), pch = 16, col = color, xlab = " ", ylab = "absolute Residuen", main = "Residuen-Plot fuer M-Schaetzer") abline(h = 0, lty = 2) for(i in 1:20){ lines(c(i,i), c(0,abs(modelR$residuals[i]))) } # (g) # Hinweis: Entfernen von Ausreißern ist für explorative Analysen interessant. Falls Qualitätsmanagement betrieben wird, # muss dies unbedingt detailliert ausdiskutiert werden! data_new <- data[-c(4,8),] modelKQ_new <- lm(y ~ x1+x2+x3+x4, data_new) summary(modelKQ_new) round(abs(modelKQ_new$coefficients - beta)/beta, digits = 2) # Bemerkung: absoluter Fehler für beta4 unsinnig, da fast bei 0 ! color <- c(rep("black", 8), "green", rep("black", 5), "blue", rep("black", 2), "cyan") plot(c(1:3, 5:7, 9:20), modelKQ_new$residuals, pch = 16, col = color, xlab = " ", ylab = "Residuen", main = "Residuen-Plot für KQ-Schaetzer") abline(h = 0, lty = 2) Mres1 <- c(modelKQ_new$residuals[1:3],0,modelKQ_new$residuals[4:6], 0, modelKQ_new$residuals[7:18]) for(i in c(1:3, 5:7, 9:20)){ lines(c(i,i), c(0,Mres1[i])) } plot(c(1:3, 5:7, 9:20), abs(modelKQ_new$residuals), pch = 16, col = color, xlab = " ", ylab = "absolute Residuen", main = "Residuen-Plot für KQ-Schaetzer") abline(h = 0, lty = 2) for(i in c(1:3, 5:7, 9:20)){ lines(c(i,i), c(0,abs(Mres1[i]))) } data_new2 <- data[-4,] modelKQ_new2 <- lm(y ~ x1+x2+x3+x4, data_new2) summary(modelKQ_new2) round(abs(modelKQ_new2$coefficients - beta)/beta, digits = 2) # Bemerkung: absoluter Fehler für beta4 unsinnig, da fast bei 0 ! color <- c(rep("black", 6), "orange", "black", "green", rep("black", 5), "blue", rep("black", 2), "cyan") plot(c(1:3, 5:8, 9:20), modelKQ_new2$residuals, pch = 16, col = color, xlab = " ", ylab = "Residuen", main = "Residuen-Plot für KQ-Schaetzer") abline(h = 0, lty = 2) Mres2 <- c(modelKQ_new2$residuals[1:3],0,modelKQ_new2$residuals[4:19]) for(i in c(1:3, 5:8, 9:20)){ lines(c(i,i), c(0,Mres2[i])) } plot(c(1:3, 5:8, 9:20), abs(modelKQ_new2$residuals), pch = 16, col = color, xlab = " ", ylab = "absolute Residuen", main = "Residuen-Plot für KQ-Schaetzer") abline(h = 0) for(i in c(1:3, 5:8, 9:20)){ lines(c(i,i), c(0,abs(Mres2[i]))) }