# Use of rnorm to simulate variance and error.
# One-predictor regression

# Use some real data for x
simRegress=function()
{
	x = cecrin$dbh1
	x = pupsize$AG
	# Or invent x
	x = 80:120
	
	
	# The model
	yhat = 20 + 1.5*x
	
	yvar = 10
	ytrue = yhat + rnorm(n=length(x),mean=0,sd=yvar)
	
	yerror = 2
	yobs = ytrue + rnorm(n=length(x),mean=0,sd=yerror)
	
	plot(x,yobs)
	fit=lm(yobs~x)
}

# Multi-level regression (two groups), each with identical model

simRegress.2Level=function()
{
	ytrue1 = yhat + rnorm(n=length(x),mean=0,sd=yvar)
	yobs1  = ytrue1 + rnorm(n=length(x),mean=0,sd=yerror)
	ytrue2 = yhat + rnorm(n=length(x),mean=0,sd=yvar)
	yobs2  = ytrue2 + rnorm(n=length(x),mean=0,sd=yerror)
	
	data1=data.frame(group='a',x=x,y=yobs1)
	data2=data.frame(group='b',x=x,y=yobs2)
	data=rbind(data1,data2)
	
	plot(data$x,data$y,pch=16,col='blue')
	bsection=data$group=='b'
	points(data$x[bsection],data$y[bsection],pch=16,col='red')
	
	fit=lm(y~x+group,data=data)
	fita=lm(y[!bsection]~x[!bsection],data=data)
	fitb=lm(y[bsection]~x[bsection],data=data)
	
	# It is nearly always best to center x on zero
	ctrx=data$x-mean(data$x)
	fit=lm(y~ctrx+group,data=data)
	fita=lm(y[!bsection]~ctrx[!bsection],data=data)
	fitb=lm(y[bsection]~ctrx[bsection],data=data)
	
	# A short-cut to adjust lines back to original x values
	linea=as.matrix(parallel.line(coef(fita)[1],coef(fita)[2],100,coef(fita)[1]))
	abline(linea,col='blue')
	lineb=as.matrix(parallel.line(coef(fitb)[1],coef(fitb)[2],100,coef(fita)[1]))
	abline(lineb,col='red')
	
	
	# A prelude to lmer
	mod=lmer(y~1+ctrx+(1+ctrx|group),data=data,verbose=TRUE)   ## This does not require arm
	display(mod)         ## This requires arm
	coef(mod)
	summary(mod)@REmat   ## Does not require arm
}

# Simulate a regression where intercept is fixed but slope a random
# variable. 
addRegressGroup=function(x,yerror,inter,slope.mean,slope.err,group)
{
 newslope=rnorm(1,mean=slope.mean,sd=slope.err)
 y=inter+newslope*x+rnorm(length(x),mean=0,sd=yerror)
 return(data.frame(group,x,y))
}


simulateMultiLevelRegression=function(x=(-20:20),within.sd,inter,slope.mean,slope.err,groups=50)
{
 for(i in 1:groups)
  {
   onegroup=addRegressGroup(x=x,yerror=within.sd,inter=inter,slope.mean=slope.mean,slope.err=slope.err,group=i)
   if(i==1) result=onegroup
   else result=rbind(result,onegroup)
  }
  
 result$group=as.factor(result$group)
 return(result)
}


# Convert coefficients from lm from a one-factor model into a dataframe
convertCoef.lm=function(fit) 
{
 cf=coef(fit)
 no.factor=length(cf)-1
 result=matrix(ncol=2,nrow=no.factor)
 
 result[1,]=cf[1:2]
 result[,2]=cf[2]
 for(i in 2:no.factor) result[i,1]=cf[1]+cf[i+1]
 
 return(result)
}

sample.multilevel.regress=function()
{
	data=simulateMultiLevelRegression(x=(-5):5,within.sd=15,inter=100,slope.mean=8,slope.err=0,groups=50)
	data=simulateMultiLevelRegression(x=(-5):5,within.sd=15,inter=100,slope.mean=8,slope.err=2,groups=50)
	
	std.fit=as.data.frame(coef(lmList(y~x|group,data=data)))
	lmer.fit=lmer(y ~ 1 + x +(1+x|group),data=data)
	linear.fit=as.data.frame(convertCoef.lm(lm(y~x+group,data=data)))

	xyplot(y~x|group,data=data,panel=function(x,y){panel.xyplot(x,y); panel.lmline(x,y)})
	plot(data$x,data$y)
	points(data$x[data$group==1],data$y[data$group==1],pch=16,col='red',cex=2)
	abline(as.matrix(std.fit[1,]),col='red')
	points(data$x[data$group==18],data$y[data$group==18],pch=16,col='blue',cex=2)
	abline(as.matrix(std.fit[18,]),col='blue')
	hist(std.fit[,2],breaks=25)
	
	for(i in 1:dim(std.fit)[1]) abline(as.matrix(std.fit[i,]),col='gray')
	for(i in 1:dim(coef(mod)[[1]])[1]) abline(as.matrix(coef(mod)[[1]][i,]),col='red')
	for(i in 1:dim(linear.fit)[1]) abline(as.matrix(linear.fit[i,]),col='green')
	abline(lm(y~x,data=data))
	
	xyplot(y~x|group,data=data,
	       std.list=as.data.frame(std.fit),coef.list=as.data.frame(coef(lmer.fit)[[1]]),
	       panel=function(..., coef.list,std.list) 
	        { 
	         panel.xyplot(...)
	         panel.abline(as.numeric(coef.list[packet.number(),]))
	         panel.abline(fixef(mod),col='blue')
	         panel.abline(as.numeric(std.list[packet.number(),]),col='red') 
	        })
	
	xyplot(y~x|group,data=data,
	       lin.list=linear.fit,std.list=std.fit,coef.list=as.data.frame(coef(lmer.fit)[[1]]),
	       panel=function(..., coef.list,std.list,lin.list) 
	        { 
	         panel.xyplot(...)
	         panel.abline(as.numeric(coef.list[packet.number(),]))
	         panel.abline(as.numeric(lin.list[packet.number(),]),col='green')
	         panel.abline(as.numeric(std.list[packet.number(),]),col='red') 
	        })
}