Krusche, Stefan: Visualisierung und Analyse multivariater Daten in der gartenbaulichen Beratung - Methodik, Einsatz und Vergleich datenanalytischer Verfahren |
Umsetzung ausgewählter Methoden in Genstat Codes
Die Abschnitte, die mit einem Kommentar (Text in Anführungszeichen) in Großbuchstaben beginnen, erfordern eine Eingabe, zum Beispiel die Quelle der Datendatei. Rechenschritte, die keine weiter Aktion des Nutzers benötigen, sind durch einen Kommentar mit Kleinbuchstaben gekennzeichnet. Die Codes werden nach ihrem Auftauchen im Text in Kapitel 2 aufgelistet und zwar in der folgenden Reihenfolge:
Thema
A-3-2
"Bipolare CA"
"OPENING THE DATASET"
open 'c:\\personal\\cyclamen\\data\\gs_cy4.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"POINTER VALUES AND MAXIMUM ON THE RATING SCALE"
"insert variate identifiers to form pointer"
pointer[values=s_ges_44,s_kno_44,s_wur_44,s_gil_44,s_wel_44,s_kra_44]data
"INSERT MAXIMUM AND MINIMUM VALUE ON THE RATING SCALE"
scalar[value=9]maxval
scalar[value=1]minval
"preliminary calculations"
calc data[]=data[]-minval
calc maxval=maxval-minval
calc nvars=nvalues(data)
calc nvals=nvalues(data[1])
calc nnvalsm1=(nvals-1)*-1
calc nvarsp1=nvars+1
calc dp_nvars=(2*nvars)
variate[nvalues=nvals;values=1...nvals]unit
text[nvalues=nvals]Betrieb
ftext unit;Betrieb
"LABELS FOR POSITIVE AND NEGATIVE POLS OF THE RATING SCALE"
"LABELS FOR THE UNITS (1) AND FOR THE VARIABLES"
text[nvalues=nvarsp1;values='Betrieb','Gesamt','Knospen','Wurzeln','Vergilbung','Welke',\
'Krankheiten']ratings
"DESCRIPTION FOR THE POLS"
text[nvalues=2;values='+','-']Pole
"LABELS FOR THE VARIABLES IN THE DOUBLED MATRIX"
text[nvalues=dp_nvars;values='ges+','kno+','wur+','gilb+','welke+','krank+',\
'ges-','kno-','wur-','gilb-','welke-','krank-']varname1;extra='Variable'
"+/- LABELS FOR THE VARIABLES IN THE DOUBLED MATRIX"
text[nvalues=nvars;values='Gesamt+/-','Knospen+/-','Wurzeln+/-',\
'Vergilbung+/-','Welke+/-','Krankheiten+/-']varname2;extra='Variable'
"SELECTION CRITERIA - 1 YES, 2 NO"
scalar[value=1]double
"double or simple matrix"
if double==1
"doubling tha data matrix"
variate[nvalues=nvals]maximum[1...nvars]
equate maxval;maximum
for ind=1...nvars
calc negdata[ind]=maximum[ind]-data[ind]
endfor
pointer[values=data[],negdata[]]combdata
variate[nvalues=dp_nvars;values=1...dp_nvars]colno
matrix[r=unit;c=colno]dop_mat
equate[oldf=!((1,#nnvalsm1)#dp_nvars,-1)]combdata;dop_mat
else
calc dp_nvars=nvars
variate[nvalues=dp_nvars;values=1...dp_nvars]colno
A-3-3
matrix[r=unit;c=colno]dop_mat
equate[oldf=!((1,#nnvalsm1)#dp_nvars,-1)]data;dop_mat
endif
"correspondence analysis"
corresp dop_mat;rowscores=rows;colscores=cols;roots=roots;rowinertia=rowin;colinertia=colin
calc rowscore[1,2]=rows$[*;1,2]
calc colscore[1,2]=cols$[*;1,2]
"calculations"
variate[nvalues=dp_nvars;values=(1...nvars)2]sorter
variate[nvalues=dp_nvars]cop1,cop2
equate colscore[1];cop1
equate colscore[2];cop2
for ind=1...nvars
variate[nvalues=2]colum1[ind],colum2[ind]
subset[condition=sorter.eq.ind]cop1;colum1[ind]
subset[condition=sorter.eq.ind]cop2;colum2[ind]
endfor
"producing the graph"
frame 1;xlower=0;xupper=0.7;ylower=0.3;yupper=1
frame 2;xlower=0.7;xupper=1;ylower=0;yupper=0.9
frame 3;xlower=0.2;xupper=1;ylower=0;yupper=0.3
pen 1;labels=Betrieb;linestyle=1
pen 2...nvarsp1;m=line;linestyle=1;symbol=2;labels=Pole
axes[equal=scale]1;style=box;xtitle='first dimension';ytitle='second dimension'
dgraph[keydescription='Objects and Variables';window=1;keywindow=3;\
title='Correspondence Analysis of Bipolar Data']\
rowscore[2],colum2[];rowscore[1],colum1[];description=ratings
"for output"
calc rootsumt=sum(roots)
scalar ro12[1...2]
equate roots;ro12
calc varex1=(ro12[1]/rootsumt)*100
calc varex2=(ro12[2]/rootsumt)*100
calc nvalro=nvalues(roots)
scalar rox[1...nvalro]
equate roots;rox
"absolute contributions"
calc colin12[1...2]=colin$[*;1,2]
calc colabcon[1]=colin12[1]/ro12[1]
calc colabcon[2]=colin12[2]/ro12[2]
calc rowin12[1...2]=rowin$[*;1,2]
calc rowabcon[1]=rowin12[1]/ro12[1]
calc rowabcon[2]=rowin12[2]/ro12[2]
calc rowcon1[1...nvalro]=(rowin$[*;1...nvalro])/rox[1...nvalro]
calc colcon1[1...nvalro]=(colin$[*;1...nvalro])/rox[1...nvalro]
calc rowcon2[1...nvalro]=(rowin$[*;1...nvalro])
calc colcon2[1...nvalro]=(colin$[*;1...nvalro])
factor[nvalues=nvalro;levels=nvalro;values=1...nvalro]PrinAxis
factor[nvalues=dp_nvars;levels=dp_nvars;values=1...dp_nvars]colvars
factor[nvalues=nvals;levels=nvals;values=1...nvals]rowvars
factor[nvalues=2;levels=2;values=1,2]Prinaxis
factor[nvalues=nvars;levels=nvars;values=1...nvars]colvars2
pen 1...20;labels=*;brush=16
axes[equal=no] 1;xmarks=*;ymarks=*;xtitle='principal component';ytitle=*
table[class=PrinAxis,rowvars]row1tab
equate rowcon2;row1tab
dbarchart[title='CUSUM Diagramm of Unit Contributions';window=1;keywindow=2;\
append=yes;keydescription='Units']row1tab
A-3-4
table[class=Prinaxis,rowvars]rowabtab
equate rowabcon;rowabtab
dbarchart[title='Absolute Unit Contributions';window=1;keywindow=2;\
append=yes;keydescription='Units']rowabtab
calc Unit_CON[1,2]=rowabcon[]*100
table[class=PrinAxis,colvars]col1tab
equate colcon2;col1tab
dbarchart[title='CUSUM diagramm of variable contributions';window=1;keywindow=2;\
append=yes;\
keydescription='Variables']col1tab;description=varname1
table[class=Prinaxis,colvars]colabtab
equate colabcon;colabtab
dbarchart[window=1;keywindow=2;\
title='Absolute Variable contributions';\
append=yes;keydescription='Variables']colabtab;description=varname1
calc Var_CON[1,2]=colabcon[]*100
variate[nvalues=dp_nvars]varcon1a,varcon2a
equate colabcon[1];varcon1a
equate colabcon[2];varcon2a
calc varcon1b=circulate(varcon1a;nvars)
calc varcon2b=circulate(varcon2a;nvars)
variate[nvalues=dp_nvars;values=1...dp_nvars]colcase
subset[condition=colcase<=nvars]varcon1a,varcon1b,varcon2a,varcon2b;vc1a,vc1b,vc2a,vc2b
calc vc1=vc1a+vc1b
calc vc2=vc2a+vc2b
pointer[values=vc1,vc2]vcpoint
pen 1...dp_nvars;colour=2...dp_nvars
table[class=Prinaxis,colvars2]col2tab
equate vcpoint;col2tab
dbarchart[title='Absolute Variable Contributions (accumulated)';window=1;keywindow=2;\
append=yes;keydescription='Variables']col2tab;description=varname2
variate[nvalues=nvals]rowinvar[1...nvalro]
equate rowcon2;rowinvar
variate[nvalues=dp_nvars]colinvar[1...nvalro]
equate colcon2;colinvar
"relative contributions"
calc inj=vsum(colinvar)
calc ink=vsum(rowinvar)
calc Var_COR[1...nvalro]=colinvar[]/inj
calc Unit_COR[1...nvalro]=rowinvar[]/ink
calc Var_QLT=Var_COR[1]+Var_COR[2]
calc Unit_QLT=Unit_COR[1]+Unit_COR[2]
"measures of polarisation"
calc zbarj[1...dp_nvars]=mean(combdata[])
calc kbarj[1...dp_nvars]=zbarj[]/maxval
calc polmj[1...dp_nvars]=1/(kbarj[]*(1-kbarj[]))
calc kij[1...dp_nvars]=combdata[]/maxval
calc m1kij[1...dp_nvars]=1-kij[]
calc inter1[1...dp_nvars]=kij[]*m1kij[]
calc inter2[1...dp_nvars]=sum(inter1[])
calc polei[1...dp_nvars]=nvals/inter2[]
calc polr[1...dp_nvars]=polmj[]/polei[]
variate[nvalues=nvars]polar1;extra='average'
variate[nvalues=nvars]polar2;extra='individual'
variate[nvalues=nvars]polar3;extra='relative'
equate polmj;polar1
equate polei;polar2
equate polr;polar3
"printed output"
print roots
print[iprint=*]'percentage of inertia explained through first principal axis',varex1
print[iprint=*]'percentage of inertia explained through second principal axis',varex2
A-3-5
print 'absolue unit contribution to the first two dimensions in %'
print[iprint=extra]Betrieb, Unit_CON[1],Unit_CON[2]
print 'absolue variable contribution to the first two dimensions in %'
print[iprint=extra]varname1, Var_CON[1],Var_CON[2]
print 'relative contributions and quality of the two-dimensional display of the units'
print[iprint=extra] Betrieb,Unit_QLT,Unit_COR[1],Unit_COR[2]
print 'relative contributions and quality of the two-dimensional display of the variables'
print[iprint=extra] varname1,Var_QLT,Var_COR[1],Var_COR[2]
print 'polarization'
print[iprint=extra]varname2,polar1,polar2,polar3
A-3-6
"Nichtlineare Biplots"
"OPENING THE DATASET"
open 'c:\\personal\\cyclamen\\data\\gs_cy4.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITION OF DATA POINTER"
pointer[values=s_ges_44,s_kno_44,s_wur_44,s_gil_44,s_wel_44,s_kra_44]data
"DISTANCE"
text[value='city']test
"preliminaries 1"
calc nvars=nvalues(data)
calc nvals=nvalues(data[1])
"UNIT NAMES"
ftext betrieb;unitname
"VARIABLE NAMES AND UNITS DESCRIPTION"
text[values='Betrieb','Gesamt','Knospen','Wurzeln','Vergilbung','Welke',\
'Krankheiten']varname
"NUMBER OF BIPLOT TRAJECTORIES TO LOOK AT"
scalar[value=6]noofvars
"STEPSIZE"
variate[nvalues=nvars;values=1,1,1,1,1,1]varsteps[1...nvars]
"MINIMUM VALUE"
variate[nvalues=nvars;values=1,1,1,1,1,1]varmins[1...nvars]
"MAXIMUM NUMBER OF MARKERS"
scalar[value=1001]maxmarks
"preliminaries 2 - finding the marker labels"
calc noofvap1=noofvars+1
calc nodesc=-1*(nvars-noofvars)
text[nvalues=noofvap1]varsname
equate[oldf=!(noofvars,#nodesc)]varname;varsname
scalar stepsize[1...nvars]
scalar minvalue[1...nvars]
equate varsteps;stepsize
equate varmins;minvalue
calc maxmam1=maxmarks-1
variate[nvalues=maxmarks;values=0...#maxmam1]case[1...nvars]
calc case[1...nvars]=case[]*stepsize[]+minvalue[]
calc vals[1...nvars]=case[]
calc varmin[1...nvars]=min(data[1...nvars])
calc varmax[1...nvars]=max(data[1...nvars])
calc range[1...nvars]=varmax[1...nvars]-varmin[1...nvars]
matrix[r=nvals;c=nvars]m
calc n_m1=-1*(nvals-1)
equate[oldf=!((1,#n_m1)#nvars,-1)]data;m
for ind=1...nvars
subset[condition=vals[ind]<=varmin[ind]]case[ind];lowcase[ind]
endfor
for ind=1...nvars
subset[condition=vals[ind]>=varmax[ind]]case[ind];upcase[ind]
endfor
calc nonvals[1...nvars]=nvalues(lowcase[1...nvars])\
+nvalues(upcase[1...nvars])-2
calc Nvals[1...nvars]=maxmam1-nonvals[1...nvars]
for ind=1...nvars
variate[nvalues=Nvals[ind]]marker[ind]
endfor
A-3-7
calc max_lowc[1...nvars]=max(lowcase[1...nvars])
calc min_upc[1...nvars]=min(upcase[1...nvars])
for ind=1...nvars
subset[condition=case[ind]>=max_lowc[ind] .and. case[ind]<=min_upc[ind]]\
vals[ind];marker[ind]
endfor
"centre data"
scal xbar[1...nvars]
calc cdata[1...nvars] = data[] - (xbar[] = mean(data[]))
calc max[1...nvars] = max(cdata[])
calc min[1...nvars] = min(cdata[])
calc rng[1...nvars] = max[] - min[]
calc maxrng = vmax(rng)
"calculating centred marker"
calc cmarker[1...nvars]=marker[]-xbar[]
"adding the means for the point of concurrency"
variate[nvalues=1;value=0]zero[1...nvars]
variate[nvalues=1]databar[1...nvars]
equate xbar;databar
for ind=1...nvars
append marker[ind],databar[ind]
append cmarker[ind],zero[ind]
sort marker[ind]
sort cmarker[ind]
endfor
"calculating number of markers"
calc nvarsm1=nvars-1
calc marval[1...nvars]=nvalues(marker[])
calc sum_mark=vsum(marval)
calc dif_mark[1...nvars]=sum_mark-marval[]
"pco"
fsim[similarity=sim]cdata[];test=#nvars(#test);range=maxrng
pco[p=roots] sim;dist=dist;lrv=l;centroid=c
calc psco[1,2]=l[1]$[*;1,2]
axes[equal=scale]1;style=box;xorigin=*;yorigin=*;xtitle='first dimension';\
ytitle='second dimension'
pen 1;labels=unitname
pen 2;linestyle=1
dgraph[keywindow=0;title='Principal Coordinates Analysis, centred data']\
psco[2];psco[1]
dmst[dimensions=2,1;title='PCO with MST, centred data']\
coordinates=l[1];similarity=sim;\
symbols=unitname
"calculating pseudo variables"
variate[nvalues=sum_mark]pseudo[1...nvars]
for ind=1...nvars
variate[nvalues=dif_mark[ind]]zeros[ind]
calc zeros[ind]=mean(cdata[ind])
equate !P(cmarker[ind],zeros[ind]);pseudo[ind]
endfor
variate[nvalues=nvars]oneval
equate marval;oneval
calc acc_one=cum(oneval)
scalar shifter[1...nvars]
equate acc_one;shifter
for ind=1...nvars
calc basis[ind]=shifter[ind]-marval[ind]
calc pseudo[ind]=circulate(pseudo[ind];basis[ind])
A-3-8
endfor
"forming new (pseudo and old data) dataset"
calc newno=sum_mark+nvals
variate[nvalues=newno]newdata[1...nvars]
for ind=1...nvars
equate !P(cdata[ind],pseudo[ind]);newdata[ind]
endfor
"calculating new distances"
fsim[similarity=newsim]newdata[];test=#nvars(#test);range=maxrng
calc newdist=1-newsim
"forming the trajectories"
symmetricmatrix[r=nvals]oldsim
calc oldsim=newsim$[!(1...nvals)]
for ind=1...nvars
calc first[ind]=nvals+shifter[ind]-marval[ind]+1
calc last[ind]=nvals+shifter[ind]
matrix[r=marval[ind];c=nvals]ps_sim[ind]
calc ps_sim[ind]=newdist$[!(first[ind]...last[ind]);!(1...nvals)]
endfor
lrv[r=nvals;c=2]l2
pco sim;dist=dist;lrv=l2;centroid=c
for ind =1...nvars
addpoints ps_sim[ind];lrv=l2;centroid=c;coordinates=trajemat[ind]
endfor
for ind=1...nvars
calc tra_vals[ind]=nvalues(marker[ind])
variate[nvalues=tra_vals[ind]]trajector[1,2][ind]
calc trajector[1,2][ind]=trajemat[ind]$[*;1,2]
endfor
"making the graph"
for ind=1...nvars
ftext marker[ind];variable[ind];decimals=0
endfor
pen 2...noofvap1;method=line;linestyle=1;symbol=0;\
labels=variable[];size=0.75;join=given
dgraph[keywindow=0;title='Interpolative Nonlinear Biplots']\
psco[2],trajector[2][];\
psco[1],trajector[1][]
"streching the trajectories"
subset[condition=marker[1].eq.xbar[1]]trajector[1][1];correct[1]
subset[condition=marker[1].eq.xbar[1]]trajector[2][1];correct[2]
scalar cct[1,2]
equate correct;cct
calc tra2[1][1...nvars]=((trajector[1][]-cct[1])*nvars)+cct[1]
calc tra2[2][1...nvars]=((trajector[2][]-cct[2])*nvars)+cct[2]
frame 1;xlower=0;xupper=0.7;ylower=0.3;yupper=1
frame 2;xlower=0.2;xupper=1;ylower=0;yupper=0.3
pen 1;method=point;symbol=1;labels=unitname;size=1
pen 2...noofvap1;method=line;symbol=0;labels=variable[];\
size=0.75;join=given;linestyle=1
dgraph[keywindow=2;keydescription='Objects and Variables';\
title='Interpolative Nonlinear Biplot (streched trajectories)']\
psco[2],tra2[2][];\
psco[1],tra2[1][];description=varsname
A-3-9
"PCA Residuen"
"NAME AND LOCATION OF DATA FILE"
open 'c:\\personal\\cyclamen\\data\\gs_cy3_3.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITION OF DATA AND PCA MODEL"
pointer[values=n23,n29,n41,k23,k29,k41,ph23,ph29,ph41,salz23,salz29,salz41]data
text[nvalues=1;value='yes']standard
scalar[value=2]model
scalar[value=0.05]alpha
calc nvals=nvalues(data[1])
calc nvars=nvalues(data)
variate[nvalues=nvals;values=1,2,3,4,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]case
"standardisation"
if standard .eqs. 'yes'
calculate variance[1...nvars]=var(data[])
calculate stddev[1...nvars]=sqrt(variance[])
for ind=1...nvars
if stddev[ind]>0
calculate zdata[ind]=(data[ind]-mean(data[ind]))/stddev[ind]
else
calculate zdata[ind]=('missing')
endif
endfor
calc ydata[1...nvars]=data[]
calc data[]=zdata[]
elsif standard .eqs. 'no'
calc ydata[1...nvars]=data[]
endif
matrix[r=nvals;c=1]rest
pcp[nroots=model;m=var]data;lrv=l;scores=score;residual=rest
calc sqrtroot=sqrt(l[2])
scalar sval[1...model]
equate sqrtroot;sval
calc yscore[1...model]=score$[*;1...model]
calc yscore[1...model]=yscore[1...model]/sval[1...model]
calc ysquare[1...model]=yscore[]**2
variate[nvalues=nvals]y[1...model]
equate ysquare;y
calc Tsquare=vsums(y)
calc Q=rest**2
"graph of residuals"
text[nvalues=nvals]unitname
ftext case;unitname
frame 1;xlower=0.1;xupper=0.9;ylower=0.1;yupper=0.9
axes 1;xtitle='Index';ytitle='Residual';xmarks=10
pen 1;labels=unitname
dgraph[title='PCA Residuals';keywindow=0]Q;case
dgraph[title='PCA Residuals';keywindow=0]Q;case
"citical values"
pcp[m=var]data;lrv=l2
scalar lj[1...nvars]
equate l2[2];lj
calc j=model+1
pointer[values=lj[j...nvars]]p1
calc theta1=vsums(p1)
calc p2[j...nvars]=p1[]**2
calc theta2=vsums(p2)
A-3-10
calc p3[j...nvars]=p1[]*p2[]
calc theta3=vsums(p3)
calc ho=1-(2*theta1*theta3)/(3*theta2**2)
calc calpha=ednormal(alpha)
if ho>0
calc calpha=calpha*-1
elsif ho<=0
calc calpha=calpha
endif
calc qa1=(calpha*sqrt(2*theta2*ho**2))/theta1
calc qa2=(theta2*ho*(ho-1))/theta1**2
calc qa3=(qa1+qa2+1)**(1/ho)
calc qa4=theta1*qa3
calc qalpha=qa4
calc malpha=1-alpha
calc n_p=nvals-model
calc fal=fed(malpha;model;n_p)
calc critT=(model*(nvals-1)/(nvals-model))*fal
"printed output"
print 'Q-values (PCA residuals)'
print unitname,Q
print[iprint=*] 'Critical value at alpha',alpha,'for Q'
print[iprint=*]qalpha
print 'T2-values'
print unitname,Tsquare
print[iprint=*]'Critical value at alpha',alpha,'for T2'
print[iprint=*]critT
A-3-11
"Velicers partielle Korrelations-Prozedur"
"OPENING THE DATA SET"
open 'c:\\personal\\cyclamen\\data\\gs_cy3_1.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITIONS"
pointer[values=salz23,salz29,salz41,ph23,ph29,ph41,n23,n29,n41,k23,k29,k41]data
text[nvalues=1;values='yes']standard
"preliminaries"
if standard .eqs. 'yes'
calc nvars=nvalues(data)
calc v[1...nvars]=var(data[1...nvars])
calc s[1...nvars]=sqrt(v[1...nvars])
calc data[1...nvars]=(data[1...nvars]-mean(data[1...nvars]))/s[1...nvars]
elsif standard .eqs. 'no'
calc nvars=nvalues(data)
endif
calc nvals=nvalues(data[1])
"calculation of the covariance matrix"
sspm[t=data[1...nvars]]ssp;ssp=sp;nunits=nu
fsspm ssp
calc df=nu-1
calc covmat=1/df*sp
"calculation of square roots of eigenvalues"
pcp[m=corr;p=score]data;lrv=l
scalar rootvar[1...nvars]
equate l['roots'];rootvar
for ind=1...nvars
if rootvar[ind].eq.0
calc svals=0
else
calc rootvar[ind]=sqrt(rootvar[ind])
endif
endfor
scalar sinval[1...nvars]
equate rootvar;sinval
"calculation of v-vectors"
matrix[r=nvars;c=1]m[1...nvars]
calc mt=t(l['vectors'])
equate mt;m
matrix[r=nvars;c=1]vvector[1...nvars]
calc vvector[]=sinval[]*m[]
"calculation of the matrices sq"
matrix[r=nvars;c=nvars]vv[1...nvars]
calc vv[]=rtproduct(vvector[];vvector[])
scalar[value=0]null
matrix[r=nvars;c=nvars]sv[0...nvars]
equate null;sv[0]
calc nvarsm1=nvars-1
calc sv[1...nvars]=sv[0...nvarsm1]+vv[1...nvars]
matrix[r=nvars;c=nvars]sq[0...nvarsm1]
calc sq[]=sv[nvars]-sv[0...nvarsm1]
"calculation of residual matrices rq"
calc nvarsneg=-1*nvars
diagonal[r=nvars]dsq[0...nvarsm1]
equate[oldf=!(1,nvarsneg)]sq[];dsq[]
A-3-12
diagonal[r=nvars]dsq_p12[0...nvarsm1]
calc dsq_p12[]=dsq[]**(-1/2)
matrix[r=nvars;c=nvars]rq[0...nvarsm1]
calc rq[]=product(dsq_p12[];sq[])
calc rq[]=product(rq[];dsq_p12[])
"calculation of fq values"
scalar riJ[0...nvarsm1]
calc riJ[]=sum(rq[]**2)-nvars
scalar fq[0...nvarsm1]
calc fq[]=riJ[]/(nvars*nvarsm1)
variate[nvalues=nvars;values=0...nvarsm1]princomp;decimals=0
variate[nvalues=nvars;values=fq[]]fqvalues
"output and presentatation of results"
frame window=1;ylower=0;yupper=1;xlower=0;xupper=1
axes window=1;xtitle='principal components';ytitle='fq values';xmarks=1
dgraph[t='Velicers fq-values';keyw=0]fqvalues;princomp
print 'Results of the partial correlation procedure (Velicer, 1976)'
print 'Principal components, fq-values'
print[clprint=*]princomp,fqvalues
A-3-13
"Kreuzvalidation"
"NAME AND LOCATION OF DATA FILE"
open 'c:\\personal\\cyclamen\\data\\gs_cy3_3.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITION OF DATA POINTER"
pointer[values=n23,k23,ph23,salz23,n29,n41,k29,k41,ph29,ph41,\
salz29,salz41]data
"STANDARDISATION"
text[nvalues=1;values='yes']standard
"preliminaries"
calc nvals=nvalues(data[1])
calc nvars=nvalues(data)
calc n=nvals*nvars
calc varn=n*nvars
calc nvarsm1=nvars-1
calc nvars2=nvarsm1*nvars
calc circ=nvals*nvarsm1
calc nvalsm1=nvals-1
calc nmnvars=n-nvars
calc n2=n*nvarsm1
"standardisation of data"
if standard .eqs. 'yes'
calculate variance[1...nvars]=var(data[])
calculate stddev[1...nvars]=sqrt(variance[])
for ind=1...nvars
if stddev[ind]>0
calculate zdata[ind]=(data[ind]-mean(data[ind]))/stddev[ind]
else
calculate zdata[ind]=('missing')
endif
endfor
calc ydata[1...nvars]=data[]
calc data[]=zdata[]
elsif standard .eqs. 'no'
calc ydata[1...nvars]=data[]
endif
"column Reduction"
"equate pointer to variate and circulate"
variate[nvalues=n]datavar
equate data;datavar
variate[nvalues=n]datvar_l[0...nvars]
calc datvar_l[0]=datavar
calc datvar_l[1...nvars]=circulate(datvar_l[0...nvarsm1];circ)
"create sorting variates"
variate[nvalues=n;values=1...n]case1
variate[nvalues=n;values=1...n]case2[0...nvars]
calc case2[1...nvars]=circulate(case2[0...nvarsm1];circ)
"create column reduced data sets"
calc limit=nvals*nvars-nvals+1
subset[condition=case1<limit]datvar_l[],case2[]
"sort column reduced data sets"
for i=1...nvars
sort[index=case2[i]]case2[i],datvar_l[i]
endfor
"create column reduced data matrices"
matrix[r=nvals;c=nvarsm1]credmat[1...nvars]
matrix[r=nvarsm1;c=nvals]tcredmat[1...nvars]
A-3-14
equate datvar_l[1...nvars];tcredmat[1...nvars]
calc credmat[]=t(tcredmat[])
delete[redefine=yes;list=exclusive]data,nvals,nvars,n,varn,nvarsm1,\
nvars2,circ,nvalsm1,nmnvars,n2,credmat,datavar
"row Reduction"
"generate sorting factors"
factor[nvalues=n;levels=nvars]f1
factor[nvalues=n;levels=nvals]f2
generate f1,f2
"create row reduced data set"
variate[nvalues=nmnvars]redvar[1...nvals]
for i=1...nvals
subset[condition=f2.ne.i]datavar
calc redvar[i]=datavar
variate[nvalues=n]datavar
equate data;datavar
endfor
"create row reduced data matrices"
matrix[r=nvars;c=nvalsm1]trredmat[1...nvals]
matrix[r=nvalsm1;c=nvars]rredmat[1...nvals]
equate redvar;trredmat
calc rredmat[]=t(trredmat[])
delete[redefine=yes;list=exclusive]data,nvals,nvars,n,varn,nvarsm1,\
nvars2,circ,nvalsm1,nmnvars,n2,credmat,rredmat
"singular value decomposition of reduced data matrices"
"SVDs for all matrices"
"hat types equal column reduced matrices"
"bar types equal row reduced matrices"
diagonal[r=nvars]Shat[1...nvars]
diagonal[r=nvals]Sbar[1...nvals]
matrix[r=nvals;c=nvarsm1]Uhat[1...nvars]
matrix[r=nvalsm1;c=nvars]Vbar[1...nvals]
svd credmat[1...nvars];l=Uhat[1...nvars];s=Shat[1...nvars]
svd rredmat[1...nvals];s=Sbar[1...nvals];r=Vbar[1...nvals]
for ind=1...nvars
variate[nvalues=nvals]uhatv[ind][1...nvarsm1]
calc tUhat[ind]=t(Uhat[ind])
equate tUhat[ind];uhatv[ind]
scalar shatsc[ind][1...nvarsm1]
equate Shat[ind];shatsc[ind]
endfor
for ind=1...nvals
variate[nvalues=nvars]vbarv[ind][1...nvars]
calc tVbar[ind]=t(Vbar[ind])
equate tVbar[ind];vbarv[ind]
scalar sbarsc[ind][1...nvars]
equate Sbar[ind];sbarsc[ind]
endfor
delete[redefine=yes;list=exclusive]data,nvals,nvars,n,varn,nvarsm1,\
nvars2,circ,nvalsm1,nmnvars,n2,uhatv,shatsc,vbarv,sbarsc
"preliminary data manipulations to calculate xijhat"
for i=1...nvars
for j=1...nvarsm1
calc part1[i][j]=uhatv[i][j]*sqrt(shatsc[i][j])
endfor
endfor
for i=1...nvals
for j=1...nvars
calc part2[i][j]=vbarv[i][j]*sqrt(sbarsc[i][j])
endfor
endfor
variate[nvalues=varn]part2var
append[newvector=part2var]part2[][]
variate[nvalues=varn;values=(#nvars2(1),#nvars(2))#nvals]sorter1
A-3-15
variate[nvalues=n2]part2va
subset[condition=sorter1.eq.1]part2var;part2va
variate[nvalues=n2;values=(1...nvars2)#nvals]sorter2
for ind=1...nvars2
subset[condition=sorter2.eq.ind]part2va;part2new[ind]
endfor
variate[nvalues=n2]part1var
append[newvector=part1var]part1[][]
variate[nvalues=n2;values=#nvals(1...nvarsm1)#nvars]sorter3
for ind=1...nvarsm1
variate[nvalues=n]part1v[ind]
subset[condition=sorter3.eq.ind]part1var;part1v[ind]
endfor
variate[nvalues=nvals]part1new[1...nvars2]
equate part1v;part1new
variate[nvalues=nvals]uhsh[1...nvarsm1][1...nvars]
variate[nvalues=nvals]sbvb[1...nvarsm1][1...nvars]
calc uhsh[][]=part1new[]
calc sbvb[][]=part2new[]
delete[redefine=yes;list=exclusive]data,nvals,nvars,n,varn,nvarsm1,\
nvars2,circ,nvalsm1,nmnvars,n2,uhsh,sbvb
"calculation of xijhat"
calc xijhat[1...nvarsm1][1...nvars]=uhsh[][]*sbvb[][]
"control"
matrix[r=nvals;c=nvars]mats
matrix[r=nvars;c=nvals]tmats
equate data;tmats
calc mats=t(tmats)
svd mats;l=U;s=S;r=V
"parity check"
calc Uo[1...nvars]=U
calc Vo[1...nvals]=V
calc Shato[1...nvars]=S
calc Sbaro[1...nvals]=S
for ind=1...nvars
variate[nvalues=nvals]Uhatvo[ind][1...nvarsm1]
calc tUo[ind]=t(Uo[ind])
equate tUo[ind];Uhatvo[ind]
scalar Shatsco[ind][1...nvarsm1]
equate Shato[ind];Shatsco[ind]
endfor
for ind=1...nvals
variate[nvalues=nvars]Vbarvo[ind][1...nvars]
calc tVo[ind]=t(Vo[ind])
equate tVo[ind];Vbarvo[ind]
scalar Sbarsco[ind][1...nvars]
equate Sbaro[ind];Sbarsco[ind]
endfor
for i=1...nvars
for j=1...nvarsm1
calc part1o[i][j]=Uhatvo[i][j]*sqrt(Shatsco[i][j])
endfor
endfor
for i=1...nvals
for j=1...nvars
calc part2o[i][j]=Vbarvo[i][j]*sqrt(Sbarsco[i][j])
endfor
endfor
variate[nvalues=varn]p2ovar
append[newvector=p2ovar]part2o[][]
variate[nvalues=varn;values=(#nvars2(1),#nvars(2))#nvals]sorter1
variate[nvalues=n2]p2ova
subset[condition=sorter1.eq.1]p2ovar;p2ova
variate[nvalues=n2;values=(1...nvars2)#nvals]sorter2
A-3-16
for ind=1...nvars2
subset[condition=sorter2.eq.ind]p2ova;p2onew[ind]
endfor
variate[nvalues=n2]p1ovar
append[newvector=p1ovar]part1o[][]
variate[nvalues=n2;values=#nvals(1...nvarsm1)#nvars]sorter3
for ind=1...nvarsm1
variate[nvalues=n]p1ov[ind]
subset[condition=sorter3.eq.ind]p1ovar;p1ov[ind]
endfor
variate[nvalues=nvals]p1onew[1...nvars2]
equate p1ov;p1onew
variate[nvalues=nvals]uhsho[1...nvarsm1][1...nvars]
variate[nvalues=nvals]sbvbo[1...nvarsm1][1...nvars]
calc uhsho[][]=p1onew[]
calc sbvbo[][]=p2onew[]
calc xijhato[1...nvarsm1][1...nvars]=uhsho[][]*sbvbo[][]
calc xijhatn[1...nvarsm1][1...nvars]=xijhat[][]
calc parity1[1...nvarsm1][1...nvars]=xijhato[][]*xijhat[][]
for i=1...nvarsm1
for j=1...nvars
calc minpar[i][j]=min(parity1[i][j])
if minpar[i][j].lt.0
restrict xijhatn[i][j];condition=parity1[i][j].lt.0
calc xijhatn[i][j]=xijhatn[i][j]*-1
restrict xijhatn[i][j]
endif
endfor
endfor
equate xijhatn[];xijhat[]
delete[redefine=yes;list=exclusive]data,nvals,nvars,n,varn,nvarsm1,\
nvars2,circ,nvalsm1,nmnvars,n2,xijhat
"calculation of PRESS"
for ind=1...nvarsm1
variate[nvalues=n]Mhat[ind]
equate xijhat[ind];Mhat[ind]
endfor
variate[nvalues=n]xij
equate data;xij
variate[nvalues=n;values=#n(0)]M0
scalar[value=0]DM0
calc press0=(sum(xij**2))/n
for ind=1...nvarsm1
calc xijhatM[ind]=M0+Mhat[ind]
calc M0=xijhatM[ind]
calc press[ind]=(sum((xijhatM[ind]-xij)**2))/n
calc DM[ind]=nvals+nvars-2*ind
calc DM0=DM0+DM[ind]
calc DR[ind]=nvars*nvalsm1-DM0
calc W[ind]=((press0-press[ind])/DM[ind])/(press[ind]/DR[ind])
calc press0=press[ind]
endfor
"output"
variate[nvalues=nvarsm1]Wvar
equate W;Wvar
variate[nvalues=nvarsm1;values=1...nvarsm1]pc
variate[nvalues=nvarsm1;values=#nvarsm1(1)]lims
pen 1;s=2
pen 2;s=0;m=l
axes 1;xtitle='Principal Components';ytitle='W-values';xmarks=1
dgraph[keywindow=0;title='PRESS statistics']Wvar,lims;pc,pc
axes 1;xtitle='Principal Components';ytitle='W-values';xmarks=1
dgraph[keywindow=0;title='PRESS statistics']Wvar,lims;pc,pc
print 'pc',' PRESS values W'
print[iprint=*]pc,Wvar
A-3-17
"Hauptkomponenten-Biplots mit Interpolations- und Prediktionsmarkern und CUSUM-Diagram"
"NAME AND LOCATION OF DATA FILE"
open 'c:\\phd\\cyclamen\\data\\gs_cy3_3.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITION OF DATA POINTER"
pointer[values=n23,k23,ph23,salz23,n29,k29,ph29,salz29,\
n41,k41,ph41,salz41]data
calc nvars=nvalues(data)
"UNIT NAMES"
text[values='1','2','3','4','6','7','8','9','10',\
'11','12','13','14','15','16','17','18','19','20']unitname
"VARIABLE NAMES AND UNITS DESCRIPTION"
text[values='N Woche 23','K Woche 23','pH Woche 23','Salz Woche 23',\
'N Woche 29','K Woche 29','pH Woche 29','Salz Woche 29',\
'N Woche 41','K Woche 41','pH Woche 41','Salz Woche 41',\
'Betriebe']varname
"NUMBER OF BIPLOT AXES TO LOOK AT"
scalar[value=4]noofvars
"STEPSIZE"
variate[nvalues=nvars;values=#nvars(1)]varsteps[1...nvars]
"MINIMUM VALUE"
variate[nvalues=nvars;values=#nvars(-3)]varmins[1...nvars]
"MAXIMUM NUMBER OF MARKERS"
scalar[value=1001]maxmarks
"STANDARDISATION"
text[nvalues=1;values='yes']standard
"preliminaries"
calc noofvap1=noofvars+1
calc nodesc=-1*(nvars-noofvars)
text[nvalues=noofvap1]varsname
equate[oldf=!(noofvars,#nodesc)]varname;varsname
scalar stepsize[1...nvars]
scalar minvalue[1...nvars]
equate varsteps;stepsize
equate varmins;minvalue
calc ydata[1...nvars]=data[]
"standardisation"
if standard .eqs. 'yes'
calculate variance[1...nvars]=var(data[])
calculate stddev[1...nvars]=sqrt(variance[])
for ind=1...nvars
if stddev[ind]>0
calculate zdata[ind]=(data[ind]-mean(data[ind]))/stddev[ind]
else
calculate zdata[ind]=('missing')
endif
endfor
calc data[]=zdata[]
elsif standard .eqs. 'no'
calc data[]=data[]
endif
"marker 1"
A-3-18
calc maxmam1=maxmarks-1
variate[nvalues=maxmarks;values=0...#maxmam1]case[1...nvars]
calc case[1...nvars]=case[]*stepsize[]+minvalue[]
calc vals[1...nvars]=case[]
calc nvars=nvalues(data)
calc varmin[1...nvars]=min(data[1...nvars])
calc varmax[1...nvars]=max(data[1...nvars])
calc range[1...nvars]=varmax[1...nvars]-varmin[1...nvars]
calc n=nvalues(data[1])
matrix[r=n;c=nvars]m
calc n_m1=-1*(n-1)
equate[oldf=!((1,#n_m1)#nvars,-1)]data;m
for ind=1...nvars
subset[condition=vals[ind]<=varmin[ind]]case[ind];lowcase[ind]
endfor
for ind=1...nvars
subset[condition=vals[ind]>=varmax[ind]]case[ind];upcase[ind]
endfor
calc nonvals[1...nvars]=nvalues(lowcase[1...nvars])\
+nvalues(upcase[1...nvars])-2
calc nvals[1...nvars]=maxmam1-nonvals[1...nvars]
for ind=1...nvars
variate[nvalues=nvals[ind]]marker[ind]
endfor
calc max_lowc[1...nvars]=max(lowcase[1...nvars])
calc min_upc[1...nvars]=min(upcase[1...nvars])
for ind=1...nvars
subset[condition=case[ind]>=max_lowc[ind] .and. case[ind]<=min_upc[ind]]\
vals[ind];marker[ind]
endfor
"marker 2"
pcp[m=variance]data;lrv=lrv;scores=sc
variate[nvalues=n]units[1...2]
calc units[]=sc$[*;1,2]
"marker 3"
calc meandata[1...nvars]=mean(data[])
calc diff[1...nvars]=marker[]-meandata[]
calc absdiff[1...nvars]=abs(marker[]-meandata[])
calc mindiff[1...nvars]=min(absdiff[])
calc nmarks[1...nvars]=nvalues(marker[])
for ind=1...nvars
variate[nvalues=nmarks[ind];values=1...nmarks[ind]]mult[ind]
subset[condition=absdiff[ind].eq.mindiff[ind]]mult[ind];factor[ind]
variate[nvalues=nmarks[ind];values=#nmarks[ind](factor[ind])]difffac[ind]
calc mue[ind]=stepsize[ind]*(mult[ind]-difffac[ind])
subset[condition=mue[ind].eq.0]diff[ind];lambda[ind]
endfor
variate[nvalues=nvars]vrhovar[1...nvars]
calc nnvarsm1=-1*(nvars-1)
equate[oldf=!((1,#nnvarsm1)#nvars,-1)]lrv[1];vrhovar
matrix[r=nvars;c=1]vrho1
matrix[r=nvars;c=1]vrho2
equate vrhovar[1];vrho1
equate vrhovar[2];vrho2
variate[nvalues=nvars]lambvar
equate lambda;lambvar
calc p11=lambvar*vrho1
calc p12=lambvar*vrho2
scalar part11[1...nvars]
scalar part12[1...nvars]
equate p11;part11
equate p12;part12
scalar vrhok1[1...nvars]
scalar vrhok2[1...nvars]
equate vrho1;vrhok1
equate vrho2;vrhok2
A-3-19
for ind=1...nvars
calc part21[ind]=mue[ind]*vrhok1[ind]
calc part22[ind]=mue[ind]*vrhok2[ind]
calc xaxes[ind]=part11[ind]+part21[ind]
calc yaxes[ind]=part12[ind]+part22[ind]
endfor
"marker 4"
calc nvarsp1=noofvars+1
for ind=1...nvars
ftext marker[ind];text=tm[ind]
endfor
axes[equal=scale]1;style=box;\
xtitle=‘first component;ytitle=‘second component
pen nvarsp1;labels=unitname;colour=1;m=point;symbol=1;size=1
pen 1...noofvars;m=line;linestyle=1;\
symbols=2;labels=tm[1...noofvars];size=0.75
frame 1;xlower=0;xupper=0.9;ylower=0.1;yupper=1
frame 2;xlower=0.9;xupper=1.3;ylower=0;yupper=0.9
dgraph[keywindow=2;keydescription='Objects and Units';\
title='Interpolation PCA Biplot']\
yaxes[1...noofvars],units[2];\
xaxes[1...noofvars],units[1];description=varsname
"marker 5"
calc meandata[1...nvars]=mean(data[])
calc diff[1...nvars]=marker[]-meandata[]
calc absdiff[1...nvars]=abs(marker[]-meandata[])
calc mindiff[1...nvars]=min(absdiff[])
calc nmarks[1...nvars]=nvalues(marker[])
for ind=1...nvars
variate[nvalues=nmarks[ind];values=1...nmarks[ind]]mult[ind]
subset[condition=absdiff[ind].eq.mindiff[ind]]mult[ind];factor[ind]
variate[nvalues=nmarks[ind];values=#nmarks[ind](factor[ind])]difffac[ind]
calc mue[ind]=stepsize[ind]*(mult[ind]-difffac[ind])
subset[condition=mue[ind].eq.0]diff[ind];lambda[ind]
endfor
variate[nvalues=nvars]vrhovar[1...nvars]
calc nnvarsm1=-1*(nvars-1)
equate[oldf=!((1,#nnvarsm1)#nvars,-1)]lrv[1];vrhovar
matrix[r=nvars;c=1]vrho1
matrix[r=nvars;c=1]vrho2
equate vrhovar[1];vrho1
equate vrhovar[2];vrho2
scalar vrhok1[1...nvars]
scalar vrhok2[1...nvars]
equate vrho1;vrhok1
equate vrho2;vrhok2
calc nvarsm1=nvars-1
variate[nvalues=nvars;values=1,#nvarsm1(0)]ej[1...nvars]
calc ej[2...nvars]=circulate(ej[1...nvarsm1])
matrix[r=1;c=nvars]ek[1...nvars]
equate ej;ek
matrix[r=nvars;c=2]vrho
equate[oldf=!((1,nnvarsm1)2,-1)]!P(vrho1,vrho2);vrho
scalar resultx1[1...nvars]
scalar resulty1[1...nvars]
for ind=1...nvars
calc divisor[ind]=ek[ind]*+vrho*+t(vrho)*+t(ek[ind])
calc resx1[ind]=vrhok1[ind]/divisor[ind]
calc resy1[ind]=vrhok2[ind]/divisor[ind]
endfor
equate resx1;resultx1
equate resy1;resulty1
variate[nvalues=nvars]resultx2,resulty2
equate resx1;resultx2
equate resy1;resulty2
variate[nvalues=nvars]lambvar
A-3-20
equate lambda;lambvar
calc pp11=lambvar*resultx2
calc pp12=lambvar*resulty2
scalar ppart11[1...nvars]
scalar ppart12[1...nvars]
equate pp11;ppart11
equate pp12;ppart12
for ind=1...nvars
calc ppart21[ind]=mue[ind]*resultx1[ind]
calc ppart22[ind]=mue[ind]*resulty1[ind]
calc pxaxes[ind]=ppart11[ind]+ppart21[ind]
calc pyaxes[ind]=ppart12[ind]+ppart22[ind]
endfor
"marker 6"
if standard .eqs. 'no'
calc nvarsp1=noofvars+1
for ind=1...nvars
ftext marker[ind];text=tm[ind]
endfor
axes[equal=scale]1;style=box;\
xtitle=‘first component;ytitle=‘second component
pen nvarsp1;labels=unitname;colour=1;symbols=1;size=1;m=point
pen 1...noofvars;m=line;linestyle=1;\
symbols=2;labels=tm[1...noofvars];size=1
frame 1;xlower=0;xupper=0.8;ylower=0.2;yupper=1
frame 2;xlower=0.85;xupper=1.3;ylower=0;yupper=0.9
dgraph[keywindow=2;keydescription='Objects and Variables';\
title='Prediction PCA Biplot']\
pyaxes[1...noofvars],units[2];\
pxaxes[1...noofvars],units[1];description=varsname
"marker 6 ZI"
elsif standard .eqs. 'yes'
calc nvarsp1=noofvars+1
for ind=1...nvars
ftext marker[ind];text=tm[ind]
endfor
axes[equal=scale]1;style=box;\
xtitle=‘first component;ytitle=‘second component
pen nvarsp1;labels=unitname;colour=1;symbols=1;size=1;m=point
pen 1...noofvars;m=line;linestyle=1;\
symbols=2;labels=tm[1...noofvars];size=1
frame 1;xlower=0;xupper=0.9;ylower=0.1;yupper=1
frame 2;xlower=0.9;xupper=1.3;ylower=0.5;yupper=0.9
dgraph[keywindow=2;keydescription='Objects and Variables';\
endaction=continue;title='Prediction PCA Biplot']\
pyaxes[1...noofvars],units[2];\
pxaxes[1...noofvars],units[1];description=varsname
variate[nvalues=2]ypoints
variate[nvalues=2]xpoints
scalar ypoint
scalar upoint[1,2]
dread[p=*;cursortype=4]y=ypoints;x=xpoints;ygiven=units[2];xgiven=units[1];\
ysave=unitsave[2];xsave=unitsave[1]
equate[oldf=!(-1,1)]ypoints;ypoint
equate unitsave[1];upoint[1]
equate unitsave[2];upoint[2]
for ind=1...nvars
model marker[ind]
fit [p=*]pyaxes[ind]
rkeep estimates=ps[ind]
endfor
scalar gradient[1...nvars]
scalar cons[1...nvars]
equate[oldf=!(-1,1)]ps[];gradient[]
A-3-21
equate ps[];cons[]
variate[nvalues=noofvars;values=1...noofvars]ident1
calc ident1=ident1*-1
variate[nvalues=noofvars;values=#noofvars(0.1)]ident2
text[nvalues=noofvars]names
equate varsname;names
axes[equal=no]3;style=none
frame 3;xlower=0.9;xupper=1.3;ylower=0;yupper=0.5
pen 1;m=point;symbol=2;size=1;rotation=-60;labels=names
dgraph[keywindow=0;window=3;title='Variable Selection';\
endaction=continue;screen=keep]ident1;ident2
variate[nvalues=1]yyy
variate[nvalues=1]xxx
dread[print=*;window=3]y=yyy;x=xxx;ygiven=ident1;xgiven=ident2
calculate yyy=round(yyy)
calc yyy=yyy*-1
scalar varsel
equate yyy;varsel
calc predval=(gradient[varsel]*ypoint)*stddev[varsel]+mean(ydata[varsel])
calc varno=nvalues(varname)
variate[nvalues=varno;values=1...varno]varcount
restrict varname;condition=varcount.eq.varsel
restrict unitname;\
condition=(upoint[1].eq.units[1]).and.(upoint[2].eq.units[2])
print 'the predicted value for variable'
print[iprint=*]varname
print 'and unit'
print[iprint=*]unitname
print 'is'
print[iprint=*]predval
restrict unitname
restrict varname
print unitname,ydata[varsel],data[varsel]
endif
"marker 7"
calc govars=vrhovar[1]**2+vrhovar[2]**2
text[nvalues=nvars]Variable
equate[oldf=!(#nvars,-1)]varname;Variable
scalar lrvscal[1...nvars]
equate lrv[2];lrvscal
calc gounits=((lrvscal[1]+lrvscal[2])/lrv[3])*100
print[iprint=*]'The first two principal components explain',gounits
print 'percent of the total variation in the data'
print 'The adequacy of fit of the variables in two dimensions is'
print[iprint=*]Variable,govars
variate[nvalues=nvars]ev[1...nvars]
equate[oldf=!((1,nnvarsm1)#nvars,-1)]lrv[1];ev
calc evsquare[1...nvars]=ev[]**2
scalar l[1...nvars]
equate lrv[2];l
calc cusum[1...nvars]=evsquare[]*l[]
factor[nvalues=nvars;levels=nvars;values=1...nvars]PrinComp
factor[nvalues=nvars;levels=nvars;values=1...nvars;labels=Variable]vars
pen 1...20;labels=*;brush=16
axes 1;ytitle='Eigenvalue';xtitle='Principal Components'
table[class=PrinComp,vars]cusumtab
equate cusum;cusumtab
dbarchart[title='CUSUM diagram';\
append=yes;keydescription='Variables']\
cusumtab
A-3-22
A-3-23
"Gruppenanlysemodell und Gamma-q-q-Plots"
"OPENING DATA SET"
open 'c:\\phd\\kennzahl\\data\\gs_kn5.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"DEFINITION OF DATA POINTER"
pointer[values=allgawp,spezp,lohnqp,lohnak,heizqm,\
eqm,glasqm,glasqmak,\
fkp,anvermp,\
beinkp,beinkak,beinkeqm,kapkoef,rdiffp,rentkoef]data
"DEFINITION OF FACTOR POINTER"
pointer[values=qm1_shn,fregion,fjahr]fdata
"USE OF COVARIANCE SSPM OR CORRELATION MATRIX"
"use corr; variance and ssp have to be checked"
text[nvalues=1;value='corr']usemat
"NORMAL OR ROBUST PCA ESTIMATORS"
"use robust or normal"
text[nvalues=1;value='robust']estim
"_______________________________________________"
"preliminaries"
calc nvars=nvalues(data)
calc nvals=nvalues(data[1])
calc fnvars=nvalues(fdata)
variate[nvalues=nvals;values=1...nvals]case
facproduct fdata;product=comb
calc lcomb=nlevels(comb)
"_______________________________________________"
"forming and storing subsets"
for ind=1...lcomb
subset[condition=comb.eq.ind]\
data[1...nvars],case;subset[1...nvars][ind],subcase[ind]
pointer[values=subset[][ind]]subspace[ind]
endfor
"_______________________________________________"
"performing pca on normal estimates"
if estim .eqs. 'normal'
for ind=1...lcomb
pcp[p=roots;m=#usemat]subspace[ind];lrv=l[ind]
lrvscree l[ind]
endfor
"performing pca on robust estimates"
elsif estim .eqs. 'robust'
for ind =1...lcomb
robsspm[p=outlier]subspace[ind];sspm=subsspm[ind]
calc subvars=nvalues(subspace[ind][1])
variate[nvalues=subvars;values=1...subvars]case_no
print case_no,subcase[ind]
pcp[p=roots;m=#usemat]subsspm[ind];lrv=l[ind]
lrvscree l[ind]
endfor
endif
A-3-24
"storing data"
open 'c:\\phd\\kennzahl\\data\\gs_subs1.glg';c=5;f=b
store[c=5;l=i;m=replace]l,subspace,subcase,\
nvars,nvals,fnvars,comb,lcomb,estim,usemat
close 5;f=b
"deletion"
delete[redefine=yes;l=e]subspace,subcase,subsspm,\
nvars,nvals,fnvars,comb,lcomb,estim,usemat
"_______________________________________________"
"groupwise analysis"
for ind =1...lcomb
print '***** Analysis of Subgroup',ind, '*****' ;decimals=0
"preliminaries"
calc subvals[ind]=nvalues(subspace[ind][1])
calc subvars[ind]=nvalues(subspace[ind])
variate[nvalues=subvals[ind]]case_no[ind]
equate subcase[ind];case_no[ind]
"NAMING VARIABLES"
calc allgawp[ind]=subset[1][ind]
calc spezp[ind]=subset[2][ind]
calc lohnqp[ind]=subset[3][ind]
calc lohnak[ind]=subset[4][ind]
calc heizqm[ind]=subset[5][ind]
calc eqm[ind]=subset[6][ind]
calc glasqm[ind]=subset[7][ind]
calc glasqmak[ind]=subset[8][ind]
calc fkp[ind]=subset[9][ind]
calc anvermp[ind]=subset[10][ind]
calc beinkp[ind]=subset[11][ind]
calc beinkak[ind]=subset[12][ind]
calc beinkeqm[ind]=subset[13][ind]
calc kapkoef[ind]=subset[14][ind]
calc rdiffp[ind]=subset[15][ind]
calc rentkoef[ind]=subset[16][ind]
"REDEFINITION OF DATA POINTER"
pointer[values=allgawp[ind],spezp[ind],lohnqp[ind],lohnak[ind],heizqm[ind],\
eqm[ind],glasqm[ind],glasqmak[ind],\
fkp[ind],anvermp[ind],\
beinkp[ind],beinkak[ind],beinkeqm[ind],kapkoef[ind],rdiffp[ind],rentkoef[ind]]data[ind]
"VARIABLE NAMES AND UNITS DESCRIPTION"
text[values='allgawp','spezp','lohnqp','lohnak','heizqm',\
'eqm','glasqm','glasqmak',\
'fkp','anvermp',\
'beinkp','beinkak','beinkeqm','kapkoef','rdiffp','rentkoef']varname
"NUMBER OF DIMENSIONS TO LOOK AT"
scalar[value=4]noofdims
"UNIT NAMES"
ftext case_no[ind];unitname[ind]
"_______________________________________________"
"Group comparisons on robust estimates"
"standardisation"
if usemat .eqs. 'corr'
text[nvalues=1;values='yes']standard
A-3-25
else
text[nvalues=1;values='no']standard
endif
if standard .eqs. 'yes'
calculate variance[ind][1...nvars]=var(data[ind][])
calculate stddev[ind][1...nvars]=sqrt(variance[ind][])
for j=1...subvars[ind]
if stddev[ind][j]>0
calculate zdata[ind][j]=(data[ind][j]-mean(data[ind][j]))/stddev[ind][j]
else
calculate zdata[ind][j]=('missing')
endif
endfor
calc data[ind][]=zdata[ind][]
elsif standard .eqs. 'no'
calc data[ind][]=data[ind][]
endif
"robust pca"
robsspm[p=outliers]data[ind];sspm=subsspm[ind]
variate[nvalues=subvals[ind];values=1...subvals[ind]]case[ind]
print case[ind],unitname[ind]
pcp[m=#usemat;nroots=noofdims]subsspm[ind];lrv=l[ind]
print l[ind][]
endfor
"deletion
delete[redefine=yes;l=e]l,lcomb,nvars,nvals,subvars,subvals,\
noofdims,varname,fdata,flevels,estim,usemat
"calculating average vectors"
calc nvarssq=nvars*nvars
matrix[r=nvars;c=nvars;values=#nvarssq(0)]H
for ind=1...lcomb
calc LtLt[ind]=rtproduct(l[ind][1];l[ind][1])
calc H=H+LtLt[ind]
endfor
variate[nvalues=nvars]varH[1...nvars]
equate H;varH
pcp[nroots=noofdims]varH;lrv=hl
svd H;singular=s
variate[nvalues=noofdims]sscos
equate s;sscos
for k=1...noofdims
calc b[k]=hl[1]$[*;k]
endfor
"calculating angles"
for k=1...noofdims
for ind=1...lcomb
calc deltat[k][ind]=arccos(sqrt(t(b[k])*+LtLt[ind]*+b[k]))*57.29578
endfor
endfor
"preparing the printout"
for k=1...noofdims
variate[nvalues=lcomb;values=deltat[k][]]delta[k]
endfor
A-3-26
variate[nvalues=lcomb;values=1...lcomb]grp
ftext grp;group
print '*** Between-Groups Comparison of Principal Components ***'
print 'average component loadings b that minimize V'
print 'directions closest to each subspace'
print varname,b[]
print 'angles formed by each group with each direction'
print group,delta[]
print 'sum of squared cosines'
print[orientation=across]sscos
print 'the groups are defined by'
print[orientation=across;iprint=*]fdata
print 'with levels'
print[iprint=*]flevels[]
"deletion
delete[redefine=yes]
"_______________________________________________"
"gamma-q-q-plots"
"DEFINITION OF EIGENVECTOR TO BE COMPARED"
scalar[value=1]ev
"forming variates from roots for boxplots"
variate[nvalues=lcomb]evals[1...nvars]
variate[nvalues=nvars]lvar[1...lcomb]
matrix[r=lcomb;c=nvars]lmat
for ind=1...lcomb
equate l[ind][2];lvar[ind]
endfor
equate lvar;lmat
calc tlmat=t(lmat)
equate tlmat;evals
"drawing the boxplots"
variate[nvalues=nvars;values=(1...nvars)]xlab
variate[nvalues=1]xsca[1...nvars]
equate xlab;xsca
ftext xsca[];pcs[1...nvars]
boxplot[title='Boxplots of eigenvalues of subgroups for \
all principal components';\
axistitle='Eigenvalue']evals[];boxlabels=pcs[]
"finding the typical vector"
calc nsets=nvalues(l)
calc nvars=sqrt(nvalues(l[1][1]))
matrix[r=nvars;c=1]evperset[1...nsets]
for ind=1...nsets
calc evperset[ind]=l[ind][1]$[*;ev]
endfor
calc tevpers[1...nsets]=t(evperset[1...nsets])
calc evpers2[1...nsets]=evperset[]*+tevpers[]
calc E=evpers2[1]
for ind=2...nsets
calc E=E+evpers2[ind]
endfor
pcp E;lrv=lall
matrix[r=nvars;c=1]evtyp
calc evtyp=lall[1]$[*;1]
A-3-27
"calculating the dissimilarities"
for ind=1...nsets
calc dissimm[ind]=t(evtyp-evperset[ind])*+(evtyp-evperset[ind])
calc dissimp[ind]=t(evtyp+evperset[ind])*+(evtyp+evperset[ind])
variate[nvalues=2]dissimv[ind]
equate !P(dissimm[ind],dissimp[ind]);dissimv[ind]
calc dissim[ind]=min(dissimv[ind])
endfor
variate[nvalues=nsets]dist
equate dissim;dist
"fitting the gamma distribution to the distances"
distribution[dist=gamma]dist;parameters=param
scalar eta,lambda
equate param;!P(eta,lambda)
variate[nvalues=nsets;values=1...nsets]index
calc p=(index-0.375)/(nsets+0.25)
calc m=eta
calc d=1/lambda
calc gamma=edgamma(p;m;d)
sort dist,index
sort gamma
"drawing the qq gamma plot"
ftext index;groups
variate[nvalues=1;values=0]null[1,2]
variate[nvalues=1]gamax[1,2]
calc gamax[1,2]=max(gamma)
variate[nvalues=2]linepoin[1,2]
append[newvector=linepoin[1]]null[1],gamax[1]
append[newvector=linepoin[2]]null[2],gamax[2]
pen 1;labels=groups;size=0.75;symbol=2
pen 2;m=line;linestyle=1;symbol=0
axes[equal=scale] 1;xtitle='Gamma quantiles';ytitle='Ordered squared distances';\
style=box
dgraph [title='Gamma q-q plot, ';key=0]dist,linepoin[2];gamma,linepoin[1]
A-3-28
"Stabilitätzsprüfung Merkmale"
"READING DATA"
UNIT [*]
FILEREAD [PRINT=summary,groups,comments,firstline;\
name='C:/personal/CYCLAMEN/DATA/ED_CY1_3.DAT';\
MISSING='*'; SEPARATOR=' '; IMETHOD=read] FGROUPS=no
pointer[values= sub_ee , abs_n1 , groe_g10 , men_g50 , sf1_mod,\
sf2_mod , bew1_kop , bew2_kop , kre_wes , sub_and , abs_g1,\
groe_w10, men_u50 , sf1_alt , sf2_alt, bew1_fus, bew2_fus,\
kre_ost]data
ftext betrieb;Betrieb
"preliminaries"
calc nvals=nvalues(data[1])
calc nvars=nvalues(data)
calc nvarsm1=nvars-1
calc nvalsm1=nvals-1
calc nvarsp1=nvars+1
calc n=nvals*nvars
calc nmnvars=n-nvars
"DESCRIPTION OF VARIABLES"
text[nvalues=nvars;value=ee ,vm1 , fg10, mg50 , sf1_m,\
sf2_m , bw1_k , bw2_k , wes , subs , vmg1,\
fw10, mw50 , sf1_a, sf2_a, bw1_f, bw2_f,\
ost]varnames
"correspondence analysis of full data matrix"
matrix[r=nvals;c=nvars]datamat
matrix[r=nvars;c=nvals]tmat
equate data;tmat
calc datamat=t(tmat)
corresp datamat;rowscore=rows;colscore=cols;\
roots=eigenf;rowinertia=rowinf;colinertia=colinf
variate[nvalues=nvals]rowf[1,2]
variate[nvalues=nvars]colf[1,2]
calc rowf[1,2]=rows$[*;1,2]
calc colf[1,2]=cols$[*;1,2]
"row reduction"
"Generate sorting factors"
factor[nvalues=n;levels=nvars]f1
factor[nvalues=n;levels=nvals]f2
generate f1,f2
"Create row reduced data set"
variate[nvalues=nmnvars]redvar[1...nvals]
for i=1...nvals
variate[nvalues=n]datavar
equate data;datavar
subset[condition=f2.ne.i]datavar
calc redvar[i]=datavar
equate data;datavar
endfor
"Create row reduced data matrices"
matrix[r=nvars;c=nvalsm1]trredmat[1...nvals]
matrix[r=nvalsm1;c=nvars]rredmat[1...nvals]
equate redvar;trredmat
calc rredmat[]=t(trredmat[])
"correspondence analysis of reduced data matrix"
A-3-29
for ind=1...nvals
corresp rredmat[ind];rowscore=row[ind];\
colscore=col[ind];roots=eigen_r[ind];\
rowinertia=rowin_r[ind];colinertia=colin_r[ind]
calc rowvar[ind][1,2]=row[ind]$[*;1,2]
calc colvar[ind][1,2]=col[ind]$[*;1,2]
endfor
"first deletion"
delete[redefine=yes;l=exclusive]varnames,betrieb,Betrieb,nvals,\
nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,cols,col,rows,row,eigen_f,eigen_r,colin_f,colin_r,\
rowin_f,rowin_r,rowf,colf,rowvar,colvar
"results from simple ca"
calc nnvalsm1=(nvals-1)*-1
calc dp_nvars=nvars
"calcs for output"
calc rootsumt=sum(eigen_f)
scalar ro12[1...2]
equate eigen_f;ro12
calc varex1=(ro12[1]/rootsumt)*100
calc varex2=(ro12[2]/rootsumt)*100
calc nvalro=nvalues(eigen_f)
scalar rox[1...nvalro]
equate eigen_f;rox
"absolute contributions"
calc colinf12[1...2]=colin_f$[*;1,2]
calc colabcon[1]=colinf12[1]/ro12[1]
calc colabcon[2]=colinf12[2]/ro12[2]
calc rowinf12[1...2]=rowin_f$[*;1,2]
calc rowabcon[1]=rowinf12[1]/ro12[1]
calc rowabcon[2]=rowinf12[2]/ro12[2]
calc rowcon1[1...nvalro]=(rowin_f$[*;1...nvalro])/rox[1...nvalro]
calc colcon1[1...nvalro]=(colin_f$[*;1...nvalro])/rox[1...nvalro]
calc rowcon2[1...nvalro]=(rowin_f$[*;1...nvalro])
calc colcon2[1...nvalro]=(colin_f$[*;1...nvalro])
calc Unit_CON[1,2]=rowabcon[]*100
calc Var_CON[1,2]=colabcon[]*100
"relative contributions"
variate[nvalues=nvals]rowinvar[1...nvalro]
equate rowcon2;rowinvar
variate[nvalues=dp_nvars]colinvar[1...nvalro]
equate colcon2;colinvar
calc inj=vsum(colinvar)
calc ink=vsum(rowinvar)
calc Var_COR[1...nvalro]=colinvar[]/inj
calc Unit_COR[1...nvalro]=rowinvar[]/ink
calc Var_QLT=Var_COR[1]+Var_COR[2]
calc Unit_QLT=Unit_COR[1]+Unit_COR[2]
"printed output"
print 'roots of full data matrix'
print eigen_f
print[iprint=*]'percentage of inertia explained \
through first principal axis',varex1
print[iprint=*]'percentage of inertia explained \
through second principal axis',varex2
print 'absolue unit contribution to the first two dimensions in %'
A-3-30
print Betrieb,Unit_CON[1],Unit_CON[2]
print 'absolue variable contribution to the first two dimensions in %'
print varnames,Var_CON[1],Var_CON[2]
print 'relative contributions and quality \
of the two-dimensional display of the units'
print Betrieb,Unit_QLT,Unit_COR[1],Unit_COR[2]
print 'relative contributions and quality \
of the two-dimensional display of the variables'
print varnames,Var_QLT,Var_COR[1],Var_COR[2]
"second deletion"
delete[redefine=yes;l=exclusive]varnames,betrieb,Betrieb,\
nvals,nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,cols,col,rows,row,eigen_f,eigen_r,colin_f,colin_r,rowin_r
"results from simple ca for reduced data matrix"
calc nnvalsm1=(nvals-1)*-1
calc dp_nvars=nvars
"calcs for output"
for ind=1...nvals
calc rrootsum[ind]=sum(eigen_r[ind])
scalar rro12[ind][1...2]
equate eigen_r[ind];rro12[ind]
calc rvarex1[ind]=(rro12[ind][1]/rrootsum[ind])*100
calc rvarex2[ind]=(rro12[ind][2]/rrootsum[ind])*100
calc rnvalro[ind]=nvalues(eigen_r[ind])
scalar rrox[ind][1...rnvalro[ind]]
equate eigen_r[ind];rrox[ind]
"absolute contributions"
calc rcolin12[ind][1...2]=colin_r[ind]$[*;1,2]
calc rcolabco[ind][1]=rcolin12[ind][1]/rro12[ind][1]
calc rcolabco[ind][2]=rcolin12[ind][2]/rro12[ind][2]
calc rrowin12[ind][1...2]=rowin_r[ind]$[*;1,2]
calc rrowabco[ind][1]=rrowin12[ind][1]/rro12[ind][1]
calc rrowabco[ind][2]=rrowin12[ind][2]/rro12[ind][2]
calc rrowcon1[ind][1...rnvalro[ind]]=(rowin_r[ind]\
$[*;1...rnvalro[ind]])/rrox[ind][1...rnvalro[ind]]
calc rcolcon1[ind][1...rnvalro[ind]]=(colin_r[ind]\
$[*;1...rnvalro[ind]])/rrox[ind][1...rnvalro[ind]]
calc rrowcon2[ind][1...rnvalro[ind]]=(rowin_r[ind]$[*;1...rnvalro[ind]])
calc rcolcon2[ind][1...rnvalro[ind]]=(colin_r[ind]$[*;1...rnvalro[ind]])
calc U_CON[ind][1,2]=rrowabco[ind][]*100
calc V_CON[ind][1,2]=rcolabco[ind][]*100
variate[nvalues=nvalsm1]rbetrieb[ind]
subset[condition=betrieb.ne.ind]betrieb;rbetrieb[ind]
"printed output 1"
print 'roots of reduced data matrix'
print eigen_r[ind]
print[iprint=*]'percentage of inertia explained \
through first principal axis',rvarex1[ind]
print[iprint=*]'percentage of inertia explained \
through second principal axis',rvarex2[ind]
print 'absolue unit contribution to the first two dimensions in %'
print rbetrieb[ind],U_CON[ind][1],U_CON[ind][2]
print 'absolue variable contribution to the first two dimensions in %'
print varnames,V_CON[ind][1],V_CON[ind][2]
"third deletion"
A-3-31
delete[redefine=yes;l=exclusive]varnames,betrieb,Betrieb,\
dp_nvars,nnvalsm1,nvals,nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,cols,col,rows,row,eigen_r,colin_r,rowin_r,rrowcon2,rcolinva,\
rcolcon2,rrowinva,rnvalro,rbetrieb,rvarex1,rvarex2
endfor
"forth deletion"
delete[redefine=yes;l=exclusive]varnames,betrieb,Betrieb,\
dp_nvars,nnvalsm1,nvals,nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,cols,col,rows,row,rrowcon2,rcolcon2,rnvalro,rbetrieb,rvarex1,rvarex2
"relative contributions"
for ind=1...nvals
variate[nvalues=nvalsm1]rrowinva[ind][1...rnvalro[ind]]
equate rrowcon2[ind];rrowinva[ind]
variate[nvalues=dp_nvars]rcolinva[ind][1...rnvalro[ind]]
equate rcolcon2[ind];rcolinva[ind]
calc rinj[ind]=vsum(rcolinva[ind])
calc rink[ind]=vsum(rrowinva[ind])
calc V_COR[ind][1...rnvalro[ind]]=rcolinva[ind][]/rinj[ind]
calc U_COR[ind][1...rnvalro[ind]]=rrowinva[ind][]/rink[ind]
calc V_QLT[ind]=V_COR[ind][1]+V_COR[ind][2]
calc U_QLT[ind]=U_COR[ind][1]+U_COR[ind][2]
"printed output 2"
print 'relative contributions and quality of \
the two-dimensional display of the units'
print rbetrieb[ind],U_QLT[ind],U_COR[ind][1],U_COR[ind][2]
print 'relative contributions and quality of \
the two-dimensional display of the variables'
print varnames,V_QLT[ind],V_COR[ind][1],V_COR[ind][2]
"fifth deletion"
delete[redefine=yes;l=exclusive]varnames,betrieb,Betrieb,\
dp_nvars,nnvalsm1,nvals,nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,cols,col,rows,row,eigen_r,colin_r,rowin_r,rrowcon2,rcolinva,\
rcolcon2,rrowinva,rnvalro,rbetrieb,rvarex1,rvarex2
endfor
"sixth deletion"
delete[redefine=yes;l=exclusive]varnames,nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,col,cols,row,rows
"procrustes rotation"
for ind=1...nvals
rotate cols;col[ind];xout=fixout[ind];yout=fitout[ind];\
residuals=resi[ind]
variate[nvalues=nvars]con_wo[ind][1,2]
calc con_wo[ind][1,2]=fitout[ind]$[*;1,2]
calc confix[ind][1,2]=fitout[ind]$[*;1,2]
endfor
"seventh deletion"
delete[redefine=yes;l=exclusive]varnames,nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,col,row,con_wo,confix,resi
"rotated jacknifed configurations"
variate[nvalues=nvars;values=1...nvars]variable
ftext variable;Variable
axes[equal=scale] 1;style=box;xtitle='first principal axis';\
ytitle='second principal axis'
A-3-32
pen 1;labels=varnames;symbol=2;size=0.75
for ind=1...nvals
dgraph[title='Configuration of variables without unit n']\
con_wo[ind][2];con_wo[ind][1]
endfor
"dotplot of residuals"
variate[nvalues=nvars]residual[1...nvals]
equate resi;residual
for ind=1...nvals
calc resimax1[ind]=max(residual[ind])
variate[nvalues=nvals]resimax2
equate resimax1;resimax2
calc resimax=max(resimax2)
endfor
for ind=1...nvals
pen 1;labels=*
if resimax1[ind]<0.05
axes[equal=no]1;xtitle=*;\
ytitle='variables';xlower=0;xupper=0.05
else
axes[equal=no]1;xtitle=*;\
ytitle='variables';xlower=0;xupper=1
endif
dotplot[g=h;title='Residuals\
of rotated jackknifed configuration of variables']\
varnames;residual[ind]
endfor
"eigth deletion"
delete[redefine=yes;l=exclusive]varnames,nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,con_wo,confix
"preparing convex hulls"
for ind =2...nvals
append con_wo[1][1],con_wo[ind][1]
append con_wo[1][2],con_wo[ind][2]
endfor
variate[nvalues=n;values=(1...nvars)#nvals]sortlist
for ind=1...nvars
variate[nvalues=nvals]peelvar[ind][1,2]
subset[condition=sortlist.eq.ind]\
con_wo[1][1],con_wo[1][2];peelvar[ind][1,2]
endfor
"ninth deletion"
delete[redefine=yes;l=exclusive]varnames,nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,confix,peelvar
for ind=1...nvars
convexhull peelvar[ind][2];peelvar[ind][1];\
yhull=yhull[ind];xhull=xhull[ind]
endfor
"drawing convex hulls"
pen 1;symbol=2;labels=varnames;size=0.5
axes[equal=scale]1;xtitle='first principal axis';\
ytitle='second principal axis';style=box;xlower=*;xupper=*
pen 2...nvarsp1;symbol=0;method=line;\
linestyle=1;join=given;colour=4
dgraph[key=0;title='Convex hulls of jackknifed \
configurations of variables']confix[1][2],yhull[1...nvars];\
confix[1][1],xhull[1...nvars]
"Stabilitätzsprüfung Objekte"
A-3-33
"READING DATA"
UNIT [*]
FILEREAD [PRINT=summary,groups,comments,firstline;\
name='C:/personal/CYCLAMEN/DATA/ED_CY1_3.DAT';\
MISSING='*'; SEPARATOR=' '; IMETHOD=read] FGROUPS=no
pointer[values= sub_ee , abs_n1 , groe_g10 , men_g50 , sf1_mod,\
sf2_mod , bew1_kop , bew2_kop , kre_wes , sub_and , abs_g1,\
groe_w10, men_u50 , sf1_alt , sf2_alt, bew1_fus, bew2_fus,\
kre_ost]data
ftext betrieb;Betrieb
"preliminaries"
calc nvals=nvalues(data[1])
calc nvars=nvalues(data)
calc nvarsm1=nvars-1
calc nvalsm1=nvals-1
calc nvarsp1=nvars+1
calc n=nvals*nvars
calc nmnvars=n-nvars
for ind =1...nvals
variate[nvalues=nvalsm1]rbetrieb[ind]
subset[condition=betrieb.ne.ind]betrieb;rbetrieb[ind]
ftext rbetrieb[ind];rBetrieb[ind]
endfor
"DESCRIPTION OF VARIABLES"
text[nvalues=nvars;value=ee ,vm1 , fg10, mg50 , sf1_m,\
sf2_m , bw1_k , bw2_k , wes , subs , vmg1,\
fw10, mw50 , sf1_a, sf2_a, bw1_f, bw2_f,\
ost]varnames
"correspondence analysis of full data matrix"
matrix[r=nvals;c=nvars]datamat
matrix[r=nvars;c=nvals]tmat
equate data;tmat
calc datamat=t(tmat)
corresp datamat;rowscore=rows;colscore=cols;\
roots=eigen_f;rowinertia=rowin_f;colinertia=colin_f
variate[nvalues=nvals]rowf[1,2]
variate[nvalues=nvars]colf[1,2]
calc rowf[1,2]=rows$[*;1,2]
calc colf[1,2]=cols$[*;1,2]
"row reduction"
"Generate sorting factors"
factor[nvalues=n;levels=nvars]f1
factor[nvalues=n;levels=nvals]f2
generate f1,f2
"Create row reduced data set"
variate[nvalues=nmnvars]redvar[1...nvals]
for i=1...nvals
variate[nvalues=n]datavar
equate data;datavar
subset[condition=f2.ne.i]datavar
calc redvar[i]=datavar
equate data;datavar
endfor
"Create row reduced data matrices"
matrix[r=nvars;c=nvalsm1]trredmat[1...nvals]
matrix[r=nvalsm1;c=nvars]rredmat[1...nvals]
equate redvar;trredmat
calc rredmat[]=t(trredmat[])
A-3-34
"correspondence analysis of reduced data matrix"
for ind=1...nvals
corresp rredmat[ind];rowscore=row[ind];\
colscore=col[ind];roots=eigen_r[ind];\
rowinertia=rowin_r[ind];colinertia=colin_r[ind]
calc rowvar[ind][1,2]=row[ind]$[*;1,2]
calc colvar[ind][1,2]=col[ind]$[*;1,2]
endfor
"first deletion"
delete[redefine=yes;l=exclusive]rbetrieb,betrieb,Betrieb,rBetrieb,nvals,\
nvars,nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,rows,row
"sixth deletion"
delete[redefine=yes;l=exclusive]rbetrieb,betrieb,Betrieb,rBetrieb,nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,row,rows
"generating fix row matrices"
variate[nvalues=nvals]rowred[1...nvarsm1]
calc rowred[1...nvarsm1]=rows$[*;1...nvarsm1]
for ind=1...nvals
variate[nvalues=nvalsm1]rrowred[ind][1...nvarsm1]
subset[condition=betrieb.ne.ind]rowred[1...nvarsm1];\
rrowred[ind][1...nvarsm1]
endfor
calc nnvalsm2=-1*(nvals-2)
matrix[r=nvalsm1;c=nvarsm1]rrows[1...nvals]
for ind =1...nvals
equate[oldf=!((1,nnvalsm2)#nvarsm1,-1)] rrowred[ind];rrows[ind]
endfor
"seventh deletion"
delete[redefine=yes;l=exclusive]rbetrieb,Betrieb,rBetrieb,betrieb,\
nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,row,rows,rrows
"procrustes rotation"
for ind=1...nvals
rotate rrows[ind];row[ind];xout=fixout[ind];yout=fitout[ind];\
residuals=resi[ind]
variate[nvalues=nvalsm1]r_con_wo[ind][1,2]
variate[nvalues=nvalsm1]r_confix[ind][1,2]
calc r_con_wo[ind][1,2]=fitout[ind]$[*;1,2]
calc r_confix[ind][1,2]=fitout[ind]$[*;1,2]
endfor
"rotated jacknifed configurations"
axes[equal=scale] 1;style=box;xtitle='first principal axis';\
ytitle='second principal axis'
for ind=1...nvals
pen 1;labels=rBetrieb[ind];symbol=1;size=0.75
dgraph[title='Configuration of units without unit n']\
r_con_wo[ind][2];r_con_wo[ind][1]
endfor
"dotplot of residuals"
variate[nvalues=nvalsm1]residual[1...nvals]
equate resi;residual
for ind=1...nvals
calc resimax1[ind]=max(residual[ind])
variate[nvalues=nvals]resimax2
equate resimax1;resimax2
calc resimax=max(resimax2)
A-3-35
endfor
for ind=1...nvals
pen 1;labels=*
if resimax1[ind]<0.01
axes[equal=no]1;xtitle=*;\
ytitle='units';xlower=0;xupper=0.01
else
axes[equal=no]1;xtitle=*;\
ytitle='units';xlower=0;xupper=0.2
endif
dotplot[g=h;title='Residuals\
of rotated jackknifed configuration of units']\
rBetrieb[ind];residual[ind]
endfor
"eigth deletion"
delete[redefine=yes;l=exclusive]Betrieb,rbetrieb,betrieb,rBetrieb,\
nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,r_con_wo,r_confix
"preparing convex hulls"
for ind =2...nvals
append r_con_wo[1][1],r_con_wo[ind][1]
append r_con_wo[1][2],r_con_wo[ind][2]
append r_confix[1][1],r_confix[ind][1]
append r_confix[1][2],r_confix[ind][2]
endfor
calc rn=nvals*nvalsm1
variate[nvalues=rn]sortlist
equate rbetrieb;sortlist
factor[nvalues=rn;levels=!(1...nvals)] f_sort
equate sortlist;f_sort
tabulate[class=f_sort] r_confix[1][1];means=rfix[1][1]
tabulate[class=f_sort] r_confix[1][2];means=rfix[1][2]
variate[nvalues=nvals]rfixvar[1][1,2]
equate rfix[1][1];rfixvar[1][1]
equate rfix[1][2];rfixvar[1][2]
for ind=1...nvals
variate[nvalues=nvalsm1]peelvar[ind][1,2]
subset[condition=sortlist.eq.ind]\
r_con_wo[1][1],r_con_wo[1][2];peelvar[ind][1,2]
endfor
"ninth deletion"
delete[redefine=yes;l=exclusive]rbetrieb,Betrieb,betrieb,rBetrieb,\
nvals,nvars,\
nvarsm1,nvalsm1,nvarsp1,n,nmvars,\
data,peelvar,rfixvar
"calculating convex hulls"
for ind=1...nvals
convexhull peelvar[ind][2];peelvar[ind][1];\
yhull=yhull[ind];xhull=xhull[ind]
endfor
"drawing convex hulls"
calc nvalsp1=nvals+1
pen 1;symbol=1;labels=Betrieb;size=0.5
axes[equal=scale]1;xtitle='first principal axis';\
ytitle='second principal axis';style=box;xlower=*;xupper=*
pen 2...nvalsp1;symbol=0;method=line;\
linestyle=1;join=given;colour=4
dgraph[key=0;title='Convex hulls of jackknifed \
configurations of units']rfixvar[1][2],yhull[1...nvals];\
rfixvar[1][1],xhull[1...nvals]
A-3-36
"Genstat-Menus zum Start aus Access in Analyse hierarchischer Liniendiagramme"
"question how to start Genstat"
QUESTION [PREAMBLE=!t('How do you want to start Genstat');\
RESPONSE=_start; \
MODE=t; DEFAULT='s'] VALUES='s','db1','db2'; CHOICE= \
'standard', \
'database kennzahlen',\
'database cyclamen'
"Nelders proc for summary statistics"
proc[par=p] 'summ' " prints summary statistics of variates "
param 'X'
getatt[nval] X; !p(nv)
scal (len,mn,v,md,mnm,mxm)[1...nv]
calc len[]=nval(#X) & mn[]=mean(#X) & v[]=var(#X)
& md[]=med(#X) & mnm[]=min(#X) & mxm[]=max(#X)
print[sq=y;ipr=*] 'length','mean',' var.','median','min.','max.'
for j=1...nv
print[sq=y;ipr=*] (len,mn,v,md,mnm,mxm)[j]
endf
endp
if _start.eq.2
delete[redefine=yes;l=e]
"opening the database"
set[in=*]
open 'c:\\phd\\kennzahl\\fba\\data\\kennzahl.glg';c=5;f=b
retrieve[c=5;l=a]
close 5;f=b
"reading data from the clipboard"
UNIT [*]
FILEREAD [PRINT=summary,groups,comments,firstline;\
name='C:/phd/kennzahl/fba/schema/cases.txt';\
MISSING='*'; SEPARATOR=' '; IMETHOD=supply]C1; FGROUPS=no
"graphical environment"
frame 1;xlower=0;xupper=1;ylower=0;yupper=1
"using question to select additional statistics"
QUESTION [PREAMBLE=!t('Further Analysis of Groups',*, \
'Do you want to see a SCATTERPLOT MATRIX'); RESPONSE=_sourc1; \
MODE=t; DEFAULT='y'] VALUES='y','n'; CHOICE= \
'yes', \
'no'
if _sourc1.eq.1
QUESTION [PREAMBLE=!t('Further Analysis of Groups - SCATTERPLOT MATRIX',*, \
'Select the variables to look at'); RESPONSE=_scvar; \
list=yes;DECLARED=yes; TYPE=variate; PRESENT=yes]
endif
QUESTION [PREAMBLE=!t('Further Analysis of Groups',*, \
'Do you want to see a BOXPLOT'); RESPONSE=_sourc3; \
MODE=t; DEFAULT='y'] VALUES='y','n'; CHOICE= \
'yes', \
'no'
if _sourc3.eq.1
A-3-37
QUESTION [PREAMBLE=!t('Further Analysis of Groups - BOXPLOTS',*, \
'Select the variables to look at'); RESPONSE=_bvar; \
list=yes;DECLARED=yes; TYPE=variate; PRESENT=yes]
endif
QUESTION [PREAMBLE=!t('Further Analysis of Groups',*, \
'Do you want to look at some SUMMARY STATISTICS'); RESPONSE=_sourc4; \
MODE=t; DEFAULT='y'] VALUES='y','n'; CHOICE= \
'yes', \
'no'
if _sourc4.eq.1
QUESTION [PREAMBLE=!t('Further Analysis of Groups - SUMMARY STATISTICS',*, \
'Select the variables to look at'); RESPONSE=_svar; \
list=yes;DECLARED=yes; TYPE=variate; PRESENT=yes]
endif
QUESTION [PREAMBLE=!t('Further Analysis of Groups',*, \
'Do you want to see a DOTPLOT'); RESPONSE=_sourc2; \
MODE=t; DEFAULT='y'] VALUES='y','n'; CHOICE= \
'yes', \
'no'
if _sourc2.eq.1
QUESTION [PREAMBLE=!t('Further Analysis of Groups - DOTPLOTS',*, \
'Select the variables to look at'); RESPONSE=_dvar; \
list=yes;DECLARED=yes; TYPE=variate; PRESENT=yes]
endif
"subsets for making the scatterplot matrix"
if _sourc1.eq.1
restrict _scvar[];condition=case.in.C1
"making the scatterplot-matrix"
dscatter _scvar[]
restrict _scvar[]
endif
"subsets for making boxplots"
if _sourc3.eq.1
restrict _bvar[];condition=case.in.C1
"making the boxplot"
calc nvars=nvalues(_bvar)
for ind=1...nvars
boxplot [g=h;m=schematic;title='boxplot for the selected group']_bvar[ind]
endfor
restrict _bvar[]
endif
"summary statistics for the whole data set"
if _sourc4.eq.1
print 'summary statistics for the whole data set'
print _svar
summ _svar[]
"summary statistics for the selected group"
restrict _svar[];condition=case.in.C1
print 'summary statistics for the selected group'
print _svar
summ _svar[]
restrict _svar[]
A-3-38
endif
"subsets for making dotplots"
if _sourc2.eq.1
ftext case;betrieb
calc nvars=nvalues(_dvar)
for ind=1...nvars
calc data[ind]=_dvar[ind]
subset[condition=case.in.C1]data[ind];_dvar[ind]
endfor
subset[condition=case.in.C1]betrieb;betriebe
"making the dotplot"
for ind =1...nvars
calc mindvar[ind]=min(_dvar[ind])
axes 1;xlower=mindvar[ind]
dotplot [g=h;direction=ascending;title='dotplot for the selected group']\
betriebe;_dvar[ind]
endfor
axes 1;xlower=*
endif
"-----------------------------------------------------------------"
elsif _start.eq.3
delete[redefine=yes;l=e]
"opening the database"
print 'cyclamen database start-up file not yet defined'
"-------------------------------------------------------------"
elsif _start.eq.1
delete[redefine=yes;l=e]
"opening the database"
print 'standard start-up'
endif
© Die inhaltliche Zusammenstellung und Aufmachung dieser Publikation sowie die elektronische Verarbeitung sind urheberrechtlich geschützt. Jede Verwertung, die nicht ausdrücklich vom Urheberrechtsgesetz zugelassen ist, bedarf der vorherigen Zustimmung. Das gilt insbesondere für die Vervielfältigung, die Bearbeitung und Einspeicherung und Verarbeitung in elektronische Systeme.
DiML DTD Version 2.0 |
Zertifizierter Dokumentenserver der Humboldt-Universität zu Berlin |
HTML - Version erstellt am: Wed May 24 16:40:53 2000 |