Y<-dget("http://www.stat.duke.edu/~pdh10/Code/hoff_2011_ba/mlcomtrade.30.30.6.10")  

hist(Y) 

#### normal scores transform 
Z<-array(qnorm( rank(c(Y),na.last="keep")/(sum(!is.na(Y))+1)) , dim(Y) )
dimnames(Z)<-dimnames(Y)
Y<-Z
Y[is.na(Y)]<- 0 
hist(Y) 
####

#### get resids from additive model
a1<-apply(Y,1,mean,na.rm=TRUE)-mean(Y,na.rm=TRUE) 
a2<-apply(Y,2,mean,na.rm=TRUE)-mean(Y,na.rm=TRUE) 
a3<-apply(Y,3,mean,na.rm=TRUE)-mean(Y,na.rm=TRUE)
a4<-apply(Y,4,mean,na.rm=TRUE)-mean(Y,na.rm=TRUE)
####

Y<-Y- mean(Y,na.rm=TRUE) - outer(a1,outer(a2,outer(a3,a4,"+"),"+"),"+") 

# 1: Make scree plots along each dimension and pick a reduced rank r

# 2: obtain the truncated and OLS hosvd with your chosen rank and compare

# 3: analyze the OLS hosvd and comment on any patterns in the data