Krusche, Stefan: Visualisierung und Analyse multivariater Daten in der gartenbaulichen Beratung - Methodik, Einsatz und Vergleich datenanalytischer Verfahren

Anhang Teil III

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

  1. Bipolare Korrespondenzanalyse
  2. Nichtlineare Biplots
  3. Hauptkomponenten Residuen
  4. Velicers partielle Korrelations-Prozedur
  5. Kreuzvalidierung
  6. Hauptkomponenten-Biplots und CUSUM-Diagram
  7. Gruppenanalysemodell und Gamma-q-q-Plots
  8. Stabilitätsprüfung
  9. Genstat-Menus zur Ergänzung der formalen Begriffsanalyse


A-3-2

Bipolare Korrespondenzanalyse

 

"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

"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

Hauptkomponenten Residuen

"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

"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

Kreuzvalidierung

"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 und CUSUM-Diagram

"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

Gruppenanalysemodell und Gamma-q-q-Plots

"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ätsprüfung

"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 zur Ergänzung der formalen Begriffsanalyse

"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


[Titelseite] [Danksagung] [1] [2] [3] [4] [5] [Bibliographie] [Abkürzungsverzeichnis] [Lebenslauf] [Selbständigkeitserklärung] [Anhang] [Anhang] [Anhang] [Anhang] [Anhang]

© 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