simulation of trios and estimation of rho

## simulation of trios and estimation of rho

## test

rm(list=ls())
options(digits=3,digits.sec=4)

source('/Volumes/iDisk/Documents/R_functions/myfunctions2.r')

pre = '/Users/haky/Clients/'
input.dir = pre %&% 'Dolan.Eileen/miscovariance/'
output.dir = input.dir
source.dir = pre %&% 'Dolan.Eileen/miscovariance/'

source(source.dir %&% "functions.r")

setwd(input.dir)

## simulate trios
## Y = mu + X + epsi 
## epsi ~ N(0,rho)
newsim = T

if(newsim)
  {
    n = 300
    m = n/3
    rho = .6
    sig2 = 10
    sig = sqrt(sig2)

    Rho3 = matrix( c(1, rho, 0, rho, 1, rho, 0, rho, 1),3,3, byrow=T)
    Rho = kronecker(diag(m), Rho3)

    X = rnorm(n)
    Z = rnorm(n)
    C = chol(Rho)
    epsi = t(C) %*% Z 

    Y = X + sig * epsi
    M = matrix( c(rep(1,n), X), n,NCOL(X)+1, byrow=F )
    triono = rep(1:m,each=3)
    sortno = rep(1:3, m)

    data = data.frame(drug=Y, gene=M[,2], triono, sortno)
    write.table(data,file=input.dir %&% "simdata.csv",sep=",",
                col.names=T,quote=F, row.names=F)
  } else {
    data = read.table(file=input.dir %&% "simdata.csv",header=T,sep=",")
  }

Y = data$drug

n = NROW(Y)
X = data$gene
M = matrix( c(rep(1,n),X), n,NCOL(X)+1, byrow=F )
m = n/3
q = NCOL(M)

system.time( optimize(ener,c(-1/sqrt(2),1/sqrt(2)),Y, M, n,m,q))
opres = optimize(ener,c(-1/sqrt(2),1/sqrt(2)),Y, M, n,m,q)

rhohat = opres$min
sig2hat = sig2fun(rhohat, Y, M, n, m, q)
betahat = betafun(rhohat, Y, M, n, m, q)

print( c(rhohat, sqrt(sig2hat), betahat ) )

## compare with simple linear regression
fit = lm(Y~M-1)
print(summary(fit)$coef)
print(summary(fit)$sigma)

## install nlme
##install.packages("nlme")
##library(nlme)
##lmefit = lme(drug~gene, data=data, corr = corExp(1,~sortno|triono), random=~1|triono)
## nlme package doesn't allow for toeplitz correlation structure