# Analysis of Dispersion according to C. R. Rao (1965) Linear Statistical # Inference and Its Applications. John Wiley & Sons, New York, 522pp. # 1) function mlreg revised from brewScor.R subroutine which solved for 1 Y. # 2) Input of a test file from Steele & Torrie Table 15.12 # 3) Make it possible for fn mlreg to print out a particular matrix (e.g. the design matrix X'X) # 4) Extract the names of variables from the attributes of yx, the input matrix. # 5) Create a unique time based name to store files in using Sys.time fn. # 6) Store some files X'Xi, X'Y, Y'Y for further use using: # 7) Finished calculating the reduction due to a factor . # 8) I redefined the use of the ask() fn and put it in the fn deff header! # 9) I Create ANOVA Table of Results #10) Factor effects are calculated. Their significance is established using df, pf, and qf. #11) I created the option to select the ny and nx indices separately. #12) Calculation of significance by subtraction is functional. #13) Extend this script to ny=1. #14) Define ndet() to calculate determinant for a scalar or 1x1 matrix. require(MASS) # define linear regression solving routine for multiple Ys and multiple Xs # YX = Y Data appended to a design matrix X # NY identifies the number of Y columns in the data set then the rest are assumed to be X # PX tells the function to print out the X'X matrix # there are ways to temporarily store results so that it will be available to routines outside the fn # rm(list=ls()) # Clean up workspace! # reduce n indices by contiguous sequence rid reduci<-function(rid, n) { ind<-1:n; n0<-length(rid); outu<- max(rid)ind; out<-sort((1:n)*(outl+outu)); out[(n0+1):n]} # ndet function ndet<-function(X) { out<-svd(X)$d l<-length(out) de<-1 for (i in 1:l) de<- de*out[i] de } # Ask function ask<- function(text) { cat(text,": ", sep="") indata<- readLines(con=stdin(), n=1, ok = TRUE) } # End of ask fn FLAM<- function(L, N, P, Q, M) { X<-rep(0,3) if ((Q==1)+(Q==2)) { if(Q==1) { df1<- P; df2<- N-P; F<-(1-L)*df2/(L*df1)} else { df1<- 2*P; df2<- 2*(N-P-1) F<- (1-L^0.5)*(N-P-1)/(P*L^0.5)} } else {if ((P==1)+(P==2)) { if(P==1) { df1<- Q; df2<- N-Q; F<-(1-L)*df2/(L*df1)} else { df1<- 2*Q; df2<- 2*(N-Q-1) F<- (1-L^0.5)*(N-Q-1)/(Q*L^0.5)} } else { s<- ((P^2 * Q^2 -4)/(P^2+Q^2-5))^0.5 la<- (P*Q-2)/4 df1<-P*Q; df2<- M*s-2*la Las<-L^(1/s) F<- (1-Las)*df2/(Las*df1) } } X<-c(df1,df2,F) X } # end of FLAM # ---------------------------------------------------------- # names(ANOVA)<-c("chisq","sq","q","n","p","s","m","df1","df2","f_ratio","PFobs","sig") # names(ANOVA)<-c("chisq","sq","q","n","p","s","m","df1","df2","f_ratio","PFobs","sig") addANOVA<- function(I, ANOV, TRAITS, CHI, S, Q, N, P, M, DF1, DF2, F, PROB, SIG) { SQ<- S*Q if (I>1) { # attach(ANOV) ANOV<-as.matrix(ANOV) CHI<- c(ANOV[,1],CHI) SQ<-c(ANOV[,2], SQ) Q<- c(ANOV[,3],Q) N<- c(ANOV[,4],N) P<- c(ANOV[,5],P) S<- c(ANOV[,6],S) M<- c(ANOV[,7],M) DF1<-c(ANOV[,8],DF1) DF2<-c(ANOV[,9],DF2) F<-c(ANOV[,10], F) PROB<- c(ANOV[,11],PROB) SIG<- c(ANOV[,12],SIG) } if (I==1) {anames<-TRAITS} else { anames<- c(row.names(ANOV),TRAITS)} rm(ANOV) # browser() Bnova<- data.frame(CHI, SQ, Q, N, P, S, M, DF1, DF2, F, PROB, SIG) row.names(Bnova)<- anames names(Bnova)<-c("chisq","sq","q","n","p","s","m","df1","df2","f_ratio","PFobs","sig") Bnova } # end of addANOVA # End of Function deffinitions # Begin Main Program # Data input crits<-c(0.05, 0.01, 0.001) SigSym<- c("ns", "*", "**", "***") print(dir(patt=".dta")) filnam<- ask("Enter a DTA filename:") cat("Sys.time:", btime<- Sys.time(), '\n') junk<- load(filnam) cat("files loaded into the workspace: \n"); print(junk) cat("all objects in the workspace: \n"); print(ls()) cat("parameters and names: \n"); print(Bsave) ny<- ncol(Bsave) noU<-rownames(Bsave)[2:nrow(Bsave)] hisdat<-Bsave[2:nrow(Bsave),] # barplot of Solution Beta (minus u estimation). barplot(t(hisdat),beside=TRUE, xlab="FACTORS", col = rainbow(ny), names.arg=noU, main="Estimated Parameters, B", legend = names(Bsave), sub="color identifies the trait affected by a factor") if (!ny==1) if (nrow(Bsave)>2) { # Normalize factors within a trait maxB<- 0*(1:ny) for (i in 1:ny) maxB[i]<- max(abs(hisdat[,i])) Nhisdat<- matrix(maxB, nrow(hisdat), ncol(hisdat), byrow=TRUE) Nhisdat<- hisdat/Nhisdat browser() # between simple and normalized solution output. barplot(t(Nhisdat),beside=TRUE, xlab="FACTORS", col = rainbow(ny), names.arg=noU, main="Normalized Estimated Parameters, B", legend = names(Bsave), sub="color identifies the trait affected by a factor") rm(Nhisdat, maxB) } else {rm(hisdat)} # end !ny==1 then nrow(Bsave)>2 item<-1:(na<-length(as.matrix(DF))) r<-as.matrix(DF)[na] df<-data.frame(DF, item) # P<- length(DiSyy[,1,1]) # not the size of set of interest fi<- 1 while(fi>0) { cat("\nAvailable Factor categories: \n"); print(t(df)) fi<- eval(parse(text=ask("Which Factor item # do you wish to test? (0 exits)"))) fy<- 1 if(fi>0) {while(sum(fy)>0) { R0<- as.matrix(DiSyy[,,na]) R1<- R0 + as.matrix(DiSyy[,,fi]) t<- r + (q<-as.matrix(DF)[fi]) cat("\nAvailable Y variables: \n"); print(t(names(Bsave))) fy<- eval(parse(text=ask("Which set of Y variables do you wish to test? (0 exits to ADISP Table)"))) f1<- 1 if ((sum(fy)>0)^(fi>0)) { ai<-0; cat("ANOVA rows set to zero\n") ANOVA<- data.frame(matrix(0,1,12)) names(ANOVA)<-c("chisq","sq","q","n","p","s","m","df1","df2","f_ratio","PFobs","sig") while(sum(f1)>0) { R0a<-R0[fy,fy]; R1a<-R1[fy,fy]; lambdaP<- ndet(R0a)/ndet(R1a); p<-length(R0a[1,]) lambdaP<- ndet(R0a)/ndet(R1a) SetOfInterest<-data.frame(fy,subnam<-names(Bsave)[fy], 1:length(fy)) names(SetOfInterest)<-c("oldID","Name ","newID") print (t(as.vector(SetOfInterest))) f1<- eval(parse(text=ask(" ... Which subset of Ys is to be tested for additional info? (0 exits subseting)"))) if(sum(f1)>0) {ai<- ai + 1; sel<-as.real(paste(setdiff(1:length(fy),f1))) s<- length(sel); R0s<- R0a[sel,sel]; R1s<-R1a[sel,sel]; lambdaS<- ndet(R0s)/ndet(R1s) lambda<- lambdaP/lambdaS n<- t; m<- n - 0.5*(p-s+q+1); chisq<- -m*log(lambda) all<-paste(fy); subnames<- gsub(" ", "",toString(subnam[f1])) cat(subnames,"/"); # cat(all,sep=',') fout<-FLAM(lambda, n, p-s, q, m) cat(' ',chisq, (p-s)*q, q, n, p, s, fout, '\n' ) prob<- pf(fout[3],fout[1], fout[2], lower.tail=FALSE) sig<-SigSym[1+sum(pf(fout[3],fout[1], fout[2], lower.tail=FALSE) 0 } } } cat("ADISP: ANALYSIS of DISPERSION for",filnam, "\n") cat("Test of additional information for", row.names(df)[fi], "\n") cat("Base traits:\n") print (t(as.vector(SetOfInterest))) cat("Additional info in:\n") ANOVA[,1]<- round(as.real(as.matrix(ANOVA[,1])),4) ANOVA[,10]<- round(as.real(as.matrix(ANOVA[,10])),4) ANOVA[,11]<- round(as.real(as.matrix(ANOVA[,11])),5) print(ANOVA) } # end of if sum fy > 0 } # end of while fi > 0 # End of program and timestamp: cat("\nSys.time:",etime<- Sys.time(), '\n') cat("Run Time:",etime-btime, 'secs.\n') rm(ai, m,n,p,q,r, R0, R0s, R0a, R1a, R1, R1s, all, ask, df, f1, fi, filnam, fy, item, junk, lambda, lambdaP, lambdaS, na, noU, SetOfInterest, chisq, FLAM, ny, reduci, s, sel, t, ANOVA, fout, sig, SigSym, crits)