# version jan 17 2010 - 3 time periods

# Graphical and informal regression analysis
# of Danish populatiom mortality rates
# 1880-84 and 2000-04 and 2005-07

raw=scan()
 0.02725     0.05213     
 0.04592     0.08235
 0.08098     0.12163 
 0.13680     0.18202 
 0.02666     0.03972
 0.04179     0.06586 
 0.06923     0.10584 
 0.11970     0.16773
 0.02359     0.03468 
 0.03934     0.05815
 0.06559     0.09622
 0.11462     0.15808
 

# Remove/add the spacing between lines 8 and 9
# if you want scan to read 24 rather than 16 #s 

r=array(raw[c( seq( 1, 7,2), seq( 2, 8,2),
               seq( 9,15,2), seq(10,16,2),
               seq(17,23,2), seq(18,24,2)  )  ], 
  c(4,2,3)) #  ages[4] sexes[2] periods[3] ;

r # 24 rates

log.r = round(log(r),4) ; 

log.r # 24 log.rates

	
# SIMPLE ADDITIVE MODEL 

# 7 parameters describe log rates in 24 cells

# Females aged 70- in 1980-84 as 'reference' cell.

# MODEL for log[rate] 
  

#  Females                         Males

#  B.0                        B.0          + B.male
#  B.0 + B.age2               B.0 + B.age2 + B.male
#  B.0 + B.age3               B.0 + B.age3 + B.male
#  B.0 + B.age4               B.0 + B.age4 + B.male
#
#  B.0          + B.20years   B.0          + B.male + B.20years
#  B.0 + B.age2 + B.20years   B.0 + B.age2 + B.male + B.20years
#  B.0 + B.age3 + B.20years   B.0 + B.age3 + B.male + B.20years
#  B.0 + B.age4 + B.20years   B.0 + B.age4 + B.male + B.20years

#
#  B.0          + B.25years   B.0          + B.male + B.25years
#  B.0 + B.age2 + B.25years   B.0 + B.age2 + B.male + B.25years
#  B.0 + B.age3 + B.25years   B.0 + B.age3 + B.male + B.25years
#  B.0 + B.age4 + B.25years   B.0 + B.age4 + B.male + B.25years


# VERY INFORMAL FITS of parameters of this model

B.hat.male   =median(log.r[,2,] - log.r[,1,]) # median of 12 estimates
B.hat.20years=median(log.r[,,2] - log.r[,,1]) # median of  8 estimates
B.hat.25years=median(log.r[,,3] - log.r[,,1]) # median of  8 estimates
B.hat.age2   =median(log.r[2,,] - log.r[1,,]) # median of  6 estimates
B.hat.age3   =median(log.r[3,,] - log.r[1,,]) # median of  6 estimates
B.hat.age4   =median(log.r[4,,] - log.r[1,,]) # median of  6 estimates

c(B.hat.male,B.hat.20years,B.hat.25years,
             B.hat.age2,B.hat.age3,B.hat.age4)

B.hat.0=median( c( log.r[1,1,1],  # arbitrary .. several other ways..
                   log.r[2,1,1]-B.hat.age2,  
                   log.r[3,1,1]-B.hat.age3,
                   log.r[4,1,1]-B.hat.age3 ) ) ; B.hat.0

# fitted values (females 1980-84, each of 4 age cateories)

fitted.log.F.1980=c(B.hat.0,
                B.hat.0+B.hat.age2,
                B.hat.0+B.hat.age3,
                B.hat.0+B.hat.age4)
# etc..

# directly  - mult. model

#ages
median( round(r[2,,]/r[1,,],2) )
median( round(r[3,,]/r[1,,],2) )
median( round(r[4,,]/r[1,,],2) )

#M:F
median( round(r[,2,]/r[,1,],2) )

#20 years later
median( round(r[,,2]/r[,,1],2) )

#25 years later
median( round(r[,,3]/r[,,1],2) )

#######################################


# homegrown median polish function

mp=function(a,i.max){
	levels=dim(a); 
	n=prod(levels)
	max.level=max(levels); dims=length(levels);
	already.removed=matrix(0,nrow=dims,ncol=max.level)
	max.abs.med=rep(1000000,dims); M.med = max(max.abs.med);
	i=0; residuals = a; #print(residuals)
	while (i < i.max & M.med > 0){
		i=i+1;
		for ( d in 1:dims) {
			med= apply(residuals,d,median);
			max.abs.med[d]=max(abs(med))
			M.med = max(max.abs.med);
			#print(med)
			if (d==1) remove = array( rep( med, n/prod(levels[-1]) ),levels)
			if (d==2) remove = array( rep( med, each=levels[1] ),levels)
			if (d==3) remove = array( rep( med, each=levels[1]*levels[2] ),levels)
			residuals = residuals - remove
			already.removed[d,1:levels[d]] = 
			  already.removed[d,1:levels[d]] + med		}
	}
	m.hat=median( already.removed[1,1:levels[1] ] )
	already.removed[1,1:levels[1] ]= 
	     already.removed[1,1:levels[1] ] - m.hat
	for (d in 1:dims) {
		if(levels[d] < max.level) 
		  already.removed[d,(levels[d]+1):max.level ] = NA
		if(d==1) fitted = already.removed[1,1:levels[1]] + m.hat;
		if(d >1) fitted = outer(fitted,already.removed[d,1:levels[d]],"+" )
		}
	result=list(effects=already.removed, 
	            observed=a,
	            fitted=fitted,
	            residuals=residuals,
	            mean=m.hat,
	            iterations=i)
	return(result)
}

f = mp(log.r,100)

f$effects
f$fitted
f$observed
f$residuals

#using corner as ref. and converting to rates

f.mult.corner1=round( exp(f$fitted - f$fitted[1,1,1]), 4)
f.mult.corner1