## Linearne hipoteze
# Napravicemo funkciju koja izvrsava testiranje linearnih hipoteza C*beta=gama
linear_hypothesis <- function(model, C, gamma) {
  # Prvo racunamo potrebne vrednosti za statistiku Q
  beta <- model$coefficients
  X <- model.matrix(model)
  cbg <- C %*% beta - gamma
  XtXi <- solve(t(X) %*% X)
  Q <- t(cbg) %*% solve(C %*% XtXi %*% t(C)) %*% cbg
  
  # Racunamo SSE
  SSE <- sum(model$residuals^2)
  n <- nrow(X)
  p <- ncol(X) - 1
  m <- nrow(C)
  
  # Test statistika je sledeca...
  Fstat <- (Q / m) / (SSE / (n-p-1))
  
  # Racunamo p vrednost na osnovu F raspodele
  pval <- 1 - pf(Fstat, m, n-p-1)
  
  # belezimo stepene slobode
  df <- c(m, n-p-1)
  
  # Ispisujemo rezultat testiranja
  cat(sprintf("F-statistic: %.4f on %d and %d DF,  p-value: %.4f\n",
                Fstat, df[1], df[2], pval))
  
  ret_val <- list("statistic" = Fstat,
                  "p.value" = pval,
                  "df" = df)
  invisible(ret_val) # vracamo listu ret_val, ali je ne prikazujemo odmah
}

# Proverimo uz jednostavan model da li su vrednosti dobre, tj iste kao u summary
# od nekog modela F statistika.
set.seed(216)
x <- runif(50, -3, 3)
z <- runif(50, -3, 3)
y <- rnorm(50)
model <- lm(y~x+z)
summary(model)
# Za test da li su svi koeficijenti nule, matrica C je
# | 0 1 0 |
# | 0 0 1 |
# a gamma je kolona vektor
# |0|
# |0|
C <- matrix(c(0,0, 1,0, 0,1), nrow=2) # matrice se popunjavaju po kolonama
gamma <- c(0,0) # u R-u su vektori kolone
linear_hypothesis(model, C, gamma)
# Iste vrednosti dobijamo kao poslednji red u summary

# Ako testiramo znacajnost koeficijenta z, gama je 0, a C je [0 0 1]
C <- t(c(0, 0, 1)) # vektor je kolona, pa moramo transponovati!
linear_hypothesis(model, C, 0)
# p vrednost je ista kao u testu za z, a dok je test statistika jednaka kvadratu
# t statistike(t value) za z
summary(model)$coef[3, 't value']^2

# I testiranje da li su u boston skupu podataka age i indus zajedno neznacajni
# se proverava linearnom hipotezom, sto smo radili anova funkcijom prosli cas
boston <- MASS::Boston
boston_model <- lm(medv ~ ., boston)
# proveravamo da li su 4. i 8. koeficijent nule odjednom
# Ako dodamo argument byrow=TRUE, u matricu se upisuje po vrstama
C <- matrix(c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
              0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), byrow = TRUE, nrow = 2)
linear_hypothesis(boston_model, C, c(0, 0))
anova(lm(medv ~ . -age - indus, boston), boston_model)
# isti rezultati kao u anova

# U primeru gde se dva koeficijenta potiru, imali smo model
set.seed(216)
x1 <- runif(50, -3, 3)
x2 <- runif(50, -3, 3)
x3 <- x2 + rnorm(50, sd = 1)
y <- 5 + x1 + x2 - x3 + rnorm(50)
model <- lm(y ~ x1 + x2 + x3)
summary(model)
# Mogli bismo da testiramo da li se koeficijenti potiru, tj da li je b2 + b3 = 0
C <- t(c(0, 0, 1, 1))
linear_hypothesis(model, C, 0)
# p vrednost je prilicno velika, sto ide u korist nasoj hipotezi da su
# koeficijenti suprotnog znaka i iste norme

# Opisimo i neki opsti princip odredjivanja matrice C i vektora gamma. Recimo da
# testiramo hipotezu H0: b1 = 1, b2 = 1, b3 = -b2 (tj. b2 + b3 = 0). Broj kolona
# matrice C mora biti broj koeficijenata vektora beta, u ovom slucaju 4. Broj
# vrsta matrice C biramo da bude jednak broju jednakosti koje imamo u hipotezi,
# dakle, sada ih je 3. Svaka vrsta matrice C treba da odredjuje jednu jednakost
# iz hipoteze, pa ako sa C_i oznacimo i-tu vrstu u C, za svaku jednakost k u H0
# treba da vazi C_k * beta = gamma_k, gde je gamma_k k-ti element vektora gamma.
# Stoga, nasa hipoteza ce matricno biti zapisana kao:
#
#             |b0|
# | 0 1 0 0 | |b1|      |1|
# | 0 0 1 0 | |b2|  =   |1|
# | 0 0 1 1 | |b3|      |0|
#
#      C   *  beta  =  gamma
C <- matrix(c(0, 1, 0, 0,
              0, 0, 1, 0,
              0, 0, 1, 1), nrow=3, byrow=TRUE)
gamma <- c(1, 1, 0)
linear_hypothesis(model, C, gamma)
# Hipoteza nije odbacena, sto smo ocekivali

# Identican rezultat dobijamo i za testiranje H0: b1 = 1, b2 = 1, b3 = -1.
# Objasniti zasto
C <- matrix(c(0, 1, 0, 0,
              0, 0, 1, 0,
              0, 0, 0, 1), nrow=3, byrow=TRUE)
gamma <- c(1, 1, -1)
linear_hypothesis(model, C, gamma)

# Za vezbu smisliti razne primere linernih hipoteza i testirati ih pomocu
# funkcije linear_hypothesis. Moguce je testirati jednakosti o bilo kakvim
# linearnim kombinacijama koeficijenata, a posto matrica C moze imati koliko god
# vrsta, moze se testirati vise jednakosti odjednom.