##Xijkl=mu+alpha i+beta j+gamma k+ alpha beta ij+alpha gamma ik+ beta gamma jk+alpha beta gamma ijk+epsilonijk

### (H0): (alpha)i = 0 for all i
### (Ha): (alpha)i != 0 for some i

### (H0): (beta)j = 0 for all j
### (Ha): (beta)j != 0 for some j

### (H0): (alpha)_ix(beta)j = 0 for all i&j
### (Ha): (alpha)_ix(beta)j != 0 for some i&j

library(GAD)
## Loading required package: matrixStats
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
observation <- c(570,1063,565,565,1080,510,583,1043,590,528,988,526,547,1026,538,521,1004,532)
P <- c(rep(1,9),rep(2,9))
tem <- rep(seq(1,3),6)
P <- as.random(P)
tem <- as.fixed(tem)
model <- aov(observation~P+tem)
GAD::gad(model)
## Analysis of Variance Table
## 
## Response: observation
##          Df Sum Sq Mean Sq  F value    Pr(>F)    
## P         1   7160    7160   16.197  0.001254 ** 
## tem       2 945342  472671 1069.257 4.924e-16 ***
## Residual 14   6189     442                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###p value is less than alpha so we reject Ho
plot(model)

# we can say var are not equal and its notmal distributed