# Dette er scriptet til oppgaven med tillegg for punktene f)-i)
# Det opprinnelige scriptet er modifisert slik at en tar vare på SSF og SSTr for hver simulering

library(nlme)
#Number of simulations
M = 2000
#Size of dataset
I = 20;J = 10
#Parameters
# punkt e)
sigma = 1;sigma.A = 0.1
# punkt f)
sigma = 1;sigma.A = 0

#Response vector and group structure
Y = rep(NA,I*J)
GR = rep(1:I,each=J)
#Vectors for results
LR = rep(NA,M)
SSE = rep(NA,M)
SSTr = rep(NA,M)
hat.sigma2 = rep(NA,M)
hat.sigma2.A = rep(NA,M)
hat.sigma2.ml = rep(NA,M)
hat.sigma2.A.ml = rep(NA,M)
#Loop for simulations
for(m in 1:M)
 {
   #Simulating data
   A = rnorm(I,0,sigma.A)
   Y = A[GR] + rnorm(I*J,0,sigma)
   #fitting full model
   fit = lme(Y~1,random=~1|GR,method="ML")
   #Fitting model without random effects
   fit2 = gls(Y~1,method="ML")
   #LR
   LR[m] = as.numeric(logLik(fit2)-logLik(fit))
   #ML estimates
   hat.sigma2.A.ml[m] = as.numeric(VarCorr(fit)[1,1])
   hat.sigma2.ml[m] = as.numeric(VarCorr(fit)[2,1])
   #Calculation of unbiased estimators
   Ybar1 = as.vector(by(Y,GR,mean))
   SSE[m] = sum((Y-rep(Ybar1,each=J))^2)
   SSTr[m] = J*sum((Ybar1-mean(Ybar1))^2)
   hat.sigma2[m] = SSE[m]/(I*(J-1))
   hat.sigma2.A[m] = SSTr[m]/(J*(I-1)) - hat.sigma2[m]/J
}


# Punkt e)
# Some plotting
par(mfrow=c(1,2))
plot(hat.sigma2,hat.sigma2.ml)
abline(c(0,1))
plot(hat.sigma2.A,hat.sigma2.A.ml)
abline(c(0,1))


# Punkt f)
# Finner antall negative verdier av hat.sigma2.A
sum(hat.sigma2.A<0)

# Skriver ut -2*LR for de simuleringene der hat.sigma2.A er negativ:
round(-2*LR[hat.sigma2.A<0],4)

# Punkt g)
# Histogram av -2*LR for de simuleringene der hat.sigma2.A er positiv:
par(mfrow=c(1,1))
hist(round(-2*LR[hat.sigma2.A>0],4),freq=FALSE,breaks=seq(0,10,0.5), xlab="",main="")

# Tegner inn tettheten til kji-kvadrat fordelingen med 1 frihets grad:
x=seq(0,10,0.1)
lines(x,dchisq(x,1))


# Punkt h)
# Finner antall forkastninger med LR testen 
sum(-2*LR>0.5*qchisq(0.95,1))


# Punkt i)
# Finner antall forkastninger med F-testen 
Fobs=(SSTr/(I-1))/(SSE/(I*(J-1)))
sum(Fobs>qf(0.95,I-1,I*(J-1)))