# # boulder workshop march 2014 # # R script # annotation: clear the memory rm(list=ls(all=TRUE)) # load OpenMx library(OpenMx) # read the data datf=read.table('rdataf') # read the female data datm=read.table('rdatam') # read the male data datb5=rbind(datf,datm) # into 1 data matrix datb5[,1]=datb5[,1]-1 # recode sex from 1->0(m) and 2->1(f) # assign variable names colnames(datb5)=c('sex', 'n1', 'n2', 'n3', 'n4', 'n5', 'n6', 'e1', 'e2', 'e3', 'e4', 'e5', 'e6', 'o1', 'o2', 'o3', 'o4', 'o5', 'o6', 'a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'c1', 'c2', 'c3', 'c4', 'c5', 'c6') # select the variables of interest isel=c(1:13) # selection of variables sex, n1-n6, e1-e6 datb2=datb5[,isel] sex=datb5[,1] # sex in a separate vector for convenience table(sex) Ss0=cov(datb2[sex==0,2:13]) # calculate the covariance matrix in males Ss1=cov(datb2[sex==1,2:13]) # calculate the covariance matrix in females Ms0=apply(datb2[sex==0,2:13],2,mean) # males means Ms1=apply(datb2[sex==1,2:13],2,mean) # females means # mi=min(Ms0,Ms1) ma=max(Ms0,Ms1) plot(1:12,Ms0,type='b',col=1,lwd=3,ylim=c(mi,ma)) lines(1:12,Ms1,type='b',col=2,lwd=3,ylim=c(mi,ma)) # # test for effect of sex on mean # alpha=.05/12; Stest=matrix(0,1,12) for (j in 2:13) { tmp=lm(datb2[,(j)]~(sex)) pval=summary(tmp)$coefficients[2,4] Stest[1,j-1]=(pval