" set working directory to current dirctory " %cd 'd:/sue/current/courses/ltcc/exercises' " load data " import 'seedsize.xls' " anova analysis " treatments line*raceme blocks block/pot/raceme anova [fpr=yes] seedsize; res=ares aplot " reml analysis: full model " vcomp [fixed=raceme] random=line+line.raceme+block/pot/raceme; \ con=none,pos,none,none,pos reml [p=#,mon] seedsize; res=rres " check residuals " vplot " look at line means " vdis [p=mean; pt=line] " save deviance " vkeep [deviance=da] " compare residuals from two analyses " dgraph rres; ares " try model without line.raceme term " vcomp [fixed=raceme] random=line+block/pot/raceme; con=pos reml [p=#,mon] seedsize " save deviance " vkeep [deviance=d0] " look at change in deviance " prin d0-da " get asymptotic p-value for change in deviance " prin 0.5-0.5*clchi(d0-da;1) " save model to check p-value " vkeep [sigma2=s2] line+block/pot; comp=vc[1,2,3] vkeep 'Constant'+raceme; eff=tc,tr scal ec; val=#tc vari er; val=tr " simulation to check p-value for line.raceme" " generate factors for interactions " facproduct [lmethod=present] !P(block,pot); bp facproduct [lmethod=present] !P(block,pot,raceme); bpr facproduct [lmethod=present] !P(line,raceme); lr calc nl,nb,np,nr = nlev(line,block,bp,bpr) " simulations " vari [nval=seedsize] yy vari [nval=10000] dev for [ntimes=10000; index=i] " print progress " if mod(i;1000).eq.0: prin [ip=*; sq=yes] i: endif " form data " calc el,eb,ebp,ebpr = grnormal(nl,nb,np,nr;0;vc[1,2,3],s2) calc yy=#ec + er$[raceme]+ el$[line] + eb$[block] + ebp$[bp] + ebpr$[bpr] " analysis " vcomp [fix=raceme] random=line+line.raceme+block/pot/raceme; con=pos reml [p=*] yy vkeep [dev=da] vcomp [fix=raceme] random=line+block/pot/raceme; con=pos reml [p=*] yy vkeep [dev=d0] " save change in deviance " calc dev$[i]=d0-da endfor " proportion of tests with zero value " prin mean(dev.gt.0) " percentiles of empirical tests distribution " prin percentiles(4(dev);90,95,99,99.5) " model without line.raceme term " vcomp [fixed=raceme] random=line+block/pot/raceme; con=pos reml [p=#,mon] seedsize " save deviance " vkeep [deviance=da] " model without line.raceme term - now try dropping line " vcomp [fixed=raceme] random=block/pot/raceme; con=pos reml [p=#,mon] seedsize " save deviance " vkeep [deviance=d0] " look at change in deviance " prin d0-da " get asymptotic p-value for change in deviance " prin 0.5-0.5*clchi(d0-da;1) " save model to check p-value " vkeep [sigma2=s2] block/pot; comp=vc[2,3] vkeep 'Constant'+raceme; eff=tc,tr scal ec; val=#tc vari er; val=tr " simulation to check p-value for line" " simulations " vari [nval=seedsize] yy vari [nval=10000] dev for [ntimes=10000; index=i] " print progress " if mod(i;1000).eq.0: prin [ip=*; sq=yes] i: endif " form data " calc eb,ebp,ebpr = grnormal(nb,np,nr;0;vc[2,3],s2) calc yy=#ec + er$[raceme]+ eb$[block] + ebp$[bp] + ebpr$[bpr] " analysis " vcomp [fix=raceme] random=line+block/pot/raceme; con=pos reml [p=*] yy vkeep [dev=da] vcomp [fix=raceme] random=block/pot/raceme; con=pos reml [p=*] yy vkeep [dev=d0] " save change in deviance " calc dev$[i]=d0-da endfor " proportion of tests with zero value " prin mean(dev.gt.0) " percentiles of empirical tests distribution " prin percentiles(4(dev);90,95,99,99.5) " final model - compare wald test with ANOVA " vcomp [fixed=raceme] random=line+block/pot/raceme; con=pos reml [p=#,mon] seedsize