#
#
#
#
#

x1 <- rnorm(20)
x2 <- rnorm(20)
x3 <- rnorm(20)

# wahres Modell
y <- 1 + 2*x1 - 1*x2 + 1*x3 + rnorm(20)
y2 <- -1 -2*x1 + x2 - x3 + rnorm(20)

X <- cbind(rep(1,20),x1,x2,x3)

# X'X
t(X) %*% X

# (X'X)^-1
solve( t(X) %*% X)

# (X'X)^-1 X'
solve(t(X) %*% X) %*% t(X)

# (X'X)^-1 X'y
solve(t(X)%*%X) %*% t(X) %*% cbind(y,y2)

# compare to 
lm(cbind(y,y2) ~ x1 + x2 + x3)

# compare to 
coef(lm(y ~ x1 + x2 + x3))
coef(lm(y2 ~ x1 + x2 + x3))


# predict
beta <- solve(t(X)%*%X) %*% t(X) %*%cbind(y,y2)

# fitted values
yhat <- X %*% beta

fitted(lm1)
summary(fitted(lm1))
summary(yhat)

summary(lm1)


x1 <- runif(100)
x2 <- runif(100)

y1 <- 1 + 2*x1 + rnorm(100)
y2 <- 1 - 2*x2 + rnorm(100)

coef(lm(y1 ~ x1 + x2))
coef(lm(y2 ~ x2 + x1))

y3 <- 1 + 2*y2 + x2*rnorm(100)
coef(lm(y3 ~ x2 + x1))
x3 <- 1 + 2*x2

lmm <- lm(cbind(y1, y2,y3) ~ x1 + x2 )

summary(lmm)

#--------------

tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 6.3,
          6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6)
gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 9.5, 9.4,
           9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2)
opacity <- c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 1.9, 5.7,
             2.8, 4.1, 3.8, 1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9)
Y <- cbind(tear, gloss, opacity)
rate     <- gl(2,10, labels = c("Low", "High"))
additive <- gl(2, 5, length = 20, labels = c("Low", "High"))

fit <- manova(Y ~ rate * additive)

lm <- lm(Y ~ rate*additive)