library(MASS)

n <- 1000
mu <- c(10,30)
sigma <- matrix(c(6,3,3,8),2,2)
mynorm1 <- mvrnorm(n, c(10,40), sigma)
mynorm2 <- cbind(mynorm1[,1]+10, mynorm1[,2]-5)
mynorm3 <- cbind(mynorm1[,1]+20, mynorm1[,2]-10)



par(mfrow = c(1,2))

plot(mean(mynorm1[,1]),mean(mynorm1[,2]), cex = 2, pch = 16, xlim = c(0,50), ylim = c(10,50), xlab = "Merkmal X", ylab = "Merkmal Y")
points(mean(mynorm2[,1]),mean(mynorm2[,2]), col = "red", cex = 2, pch = 16)
points(mean(mynorm3[,1]),mean(mynorm3[,2]), col = "blue", cex = 2, pch = 16)
all <- rbind(c(mean(mynorm1[,1]),mean(mynorm1[,2])),
	c(mean(mynorm2[,1]),mean(mynorm2[,2])),
	c(mean(mynorm3[,1]),mean(mynorm3[,2])))
abline(lm(all[,2] ~ all[,1]), lwd = 2, lty = 2)



plot(mynorm1, xlim = c(0,50), ylim = c(10,50), 
	xlab = "Merkmal X", ylab = "Merkmal Y")
points(mynorm2, col = "red")
points(mynorm3, col = "blue")
all <- rbind(mynorm1, mynorm2, mynorm3)
abline(lm(all[,2] ~ all[,1]), lwd = 2, lty = 2)

abline(lm(mynorm1[,2] ~ mynorm1[,1]), col = "black", lwd = 2)
abline(lm(mynorm2[,2] ~ mynorm2[,1]), col = "red", lwd = 2)
abline(lm(mynorm3[,2] ~ mynorm3[,1]), col = "blue", lwd = 2)

legend(10,15,c("\326kologischer Fehlschluss","'wahrer' Zusammenhang"), lty = c(2,1))

##all <- cbind(all,c(rep(1,1000),rep(2,1000),rep(3,1000)))
##summary(lm(all[,2] ~ all[,1] ))
##summary(lm(all[,2] ~ all[,1] + (all[,3])))