# Kahden selittäjän regressiomallin estimointia ja 
# hypotesin testausta havainnollistava ohjelma.
# Mallissa on jatkuva selittäjä ja yksui kaksiluokkainen luokitteluasteikollinen selittäjä.
# Todessa mallissa on interaktio selittäjien välillä,
# ja sellaisen olemassaolo oletetaan estimoitaessa.
# Ohjelma generoi simuloidun datan oletetusta todesta mallista
# ja sovittaa lineaarista regressiomallia 
# n ensimmäiseen havaintoon, kun n = 2,...,1000
# Kustakin vaiheesta esitetään parametriestimaatit, regressiokertoimen
# testin p-arvo (H0: beta_2=0) ja luottamusvälit.
# 
# Author: Lauri Mehtätalo 7.3.2017
# Saa muokata ja levittää vapaasti kunhan koodin lähde kerrotaan.  
###############################################################################

quickly<-FALSE # change this to TRUE if you want to run the simulation fast
beta<-c(-2,0.5,0.8,0.3)
windows(width=12,height=10)
par(ask=TRUE,bty="n",mai=c(1,1,0,0))
plot(1,xlim=c(0,10),ylim=c(-6,15),type="n",xlab="x",ylab="y")
abline(beta[c(1,2)],col="blue",lwd=2)
abline(beta[c(1,2)]+beta[c(3,4)],col="red",lwd=2)
nrep<-1000
X<-cbind(1,runif(nrep,0,10),c(0,1))
X<-cbind(X,X[,2]*X[,3])
y<-X%*%beta+rnorm(nrep,mean=0,sd=2)
x1<-X[,2]
x2<-X[,3]
mod<-lm(y[1:4]~x1[1:4]*x2[1:4])
ekakerta<-smallpFirstTime<-FALSE

for (i in 5:nrep) {
	
	polygon(c(0,0,12,12,0),c(8.5,15,15,8.5,8.5),col="white",border="NA")
    text(0,14.4,"Tosi malli",pos=4,font=2)
    text(0,13.6,expression(y[i]==beta[1]+beta[2]*x[i]+beta[3]*z[i]+beta[4]*x[i]*z[i]+e[i]),pos=4)
    text(0.05,12.8,expression("  "==-2+0.5*x[i]+0.8*z[i]+0.3*x[i]*z[i]+e[i]*" (tunnettu)"),pos=4)
    text(0,12,expression(var(e[i])==2^2*" kaikille i (tunnettu)" ),pos=4)
    text(0,11.2,expression(cov(e[i],e[j])==0*" kaikille i,j (tunnettu)"),pos=4)


    text(0,10,expression(H[0]*": "*beta[4]==0),pos=4,col="red")
    text(2,10,expression(H[1]*": "*beta[4]!=0),pos=4,col="red")

    text(0,9.2,expression("Kertoimen "*beta[4]*" "*95*"%"*"luottamusvali"),pos=4,col="red")
    x0<-c(0,10)

	modnew<-lm(y[1:i]~x1[1:i]*x2[1:i])
	betahat<-round(coef(modnew),1)
	emod <- bquote("  "==.(betahat[1])+.(betahat[2])*x[i]+.(betahat[3])*z[i]+.(betahat[4])*x[i]*z[i]+e[i]*" (estimoitu)")	
	if (i==2) {
       evar <-bquote(var(e[i])=="NA kaikille i.")
	   pval <-bquote(p=="NA")
	   } else {
	   evar <- bquote(var(e[i])==.(round(summary(modnew)$sigma,1))^2*" kaikille i  (estimoitu).")
	   pval <- bquote(p==.(round(summary(modnew)$coef[4,4],4)))
#	   print(anova(modnew))
       }
	   
	   if (i<4) {
		   cint <-bquote("NA,NA")
	   } else {
		   cint <- bquote(.(round(confint(modnew)[4,1],2))*", "*.(round(confint(modnew)[4,2],2)))
		   if (!ekakerta) {
			  smallpFirstTime<-summary(modnew)$coef[4,4]<0.025
		  } else {
			  smallpFirstTime<-FALSE  
		  }
		  if (!ekakerta & smallpFirstTime) ekakerta<-TRUE
	   }
	   
	text(5,14.4,paste("Estimoitu malli (OLS), n=",i),pos=4,font=2)
	text(5,13.6,expression(y[i]==beta[1]+beta[2]*x[i]+beta[3]*z[i]+beta[4]*x[i]*z[i]+e[i]),pos=4)
	text(5.05,12.8,emod,pos=4)
	text(5,12,evar,pos=4)
	text(5,11.2,expression(cov(e[i],e[j])==0*" kaikille i,j  (oletettu)"),pos=4)
	
	
	text(5,10,pval,pos=4,col="red")
	text(5,9.2,cint,pos=4,col="red")
#	xlab <- bquote(y[i]==.(betahat[1])+.(betahat[2])*x[i]^{(1)}+.(betahat[3])*x^{(2)}[i]+e[i])	
	lines(x0,coef(mod)[1]+coef(mod)[2]*x0,col="white",lwd=3)
	lines(x0,coef(mod)[1]+coef(mod)[3]+(coef(mod)[2]+coef(mod)[4])*x0,col="white",lwd=3)

	abline(beta[c(1,2)],col="blue",lwd=2)
    abline(beta[c(1,2)]+beta[c(3,4)],col="red",lwd=2)
	points(x1[i-1],y[i-1],col="white",pch=17,cex=2)
	points(x1[1:i],y[1:i],col=c("blue","red")[x2[1:i]+1])
	points(x1[i],y[i],col=c("blue","red")[x2[i]+1],pch=17,cex=2)

	lines(x0,coef(modnew)[1]+coef(modnew)[2]*x0,col="blue",lwd=2,lty="dashed")
	lines(x0,coef(modnew)[1]+coef(modnew)[3]+(coef(modnew)[2]+coef(modnew)[4])*x0,col="red",lwd=2,lty="dashed")
#	abline(modnew,col="red",lwd=2)
	mod<-modnew
	if (!quickly) {
	if (i<500) { 
		Sys.sleep(exp(-i/100)+0.1)
	} else if (!i%%10) {
		Sys.sleep(0.2)
	}
	if (i<7||i==10|i==20|i==30|i==50|i==100|i==200|i==500|smallpFirstTime) {
		text(5,-6,"Paina hiirella kuvan aluella jatkaaksesi")
		locator(1)
		polygon(c(2,2,8,8,2),c(-6.5,-5.5,-5.5,-6.5,-6.5),col="white",border="NA")
	}
}
}