## this script executes the multivariate linear Procrustes ANOVA model decomposing the shape data of Halimeda tuna segments into different sources
# first, the segments are symmetrised across the axis of their bilateral symmetry
# the dataset consists of two independent digitisations conducted in CW and CCW directions that were performed by two independent persons;thus, there are four repetitions of the original data to account for two possible sources of measurement error 
# shape variation is explained by the factors of "position", the nested spatial factors ("region", "locality" and "plant"), and their interactions

## works fine with R, ver. 4.0.2, and compatible versions of the packages uploaded below

## the script was designed for the study of Jiri Neustupa and Yvonne Nemcova entitled "..."

## uploading necessary libraries and setting the working directory for upload and saving the files
library(geomorph) ; library(shapes) ; library(MASS) ; library(pracma)
setwd("c:/R_prac/")

## this line indicates the number of repetitions in permutation tests
# thus, this needs to be indicated by a user
pocetpermutaci<-999


## now, the data of the 1st digitising person are processed

# this part will convert the original CCW coordinates along the outlines to CW and save them as a plain TXT file

pp<-readland.tps(file = "JN_CCW_IMPV.tps", readcurves = T) ; qq<-two.d.array(pp)
pocetLM<-ncol(qq) ; pocetobj<-nrow(qq) ; f<-ncol(qq) ; g<-ncol(qq)

qq2<-as.data.frame(qq)

a2<-NULL

f<-ncol(qq)
while (f>0) {
b<-qq2[,(f-1):f] ; b2<-as.matrix(b)
a2<-cbind(a2,b2)
f<-f-2 }

write.table(a2, file = "ccw_presunuto.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

# this part saves the original CW data as a plain TXT file

d<-readland.tps(file = "JN_CW_IMPV.tps", readcurves = T) ; dd<-two.d.array(d)

write.table(dd, file = "cw.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

# this part joins the CW and the (rearranged) CCW data into a single file

u<-read.table(file = "cw.txt") ; v<-read.table(file = "ccw_presunuto.txt")
JN_mat12<-rbind(u,v)

# write.table(JN_mat12, file = "JN_mat12.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

file.remove("ccw_presunuto.txt") ; file.remove("cw.txt")


## now, the data of the 2nd digitising person are processed

# this part converts the original CCW coordinates along the outlines to CW and save them as a plain TXT file

pp<-readland.tps(file = "YN_CCW_IMPV.tps", readcurves = T) ; qq<-two.d.array(pp)
pocetLM<-ncol(qq) ; pocetobj<-nrow(qq) ; f<-ncol(qq) ; g<-ncol(qq)

qq2<-as.data.frame(qq)

a2<-NULL

f<-ncol(qq)
while (f>0) {
b<-qq2[,(f-1):f] ; b2<-as.matrix(b)
a2<-cbind(a2,b2)
f<-f-2 }

write.table(a2, file = "ccw_presunuto.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

# this part saves the original CW data as a plain TXT file

d<-readland.tps(file = "YN_CW_IMPV.tps", readcurves = T) ; dd<-two.d.array(d)

write.table(dd, file = "cw.txt", append = TRUE, quote = FALSE, sep = " ",
eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

# this part joins the CW and the (rearranged) CCW data into a single file

u<-read.table(file = "cw.txt") ; v<-read.table(file = "ccw_presunuto.txt")
YN_mat12<-rbind(u,v)

# write.table(YN_mat12, file = "YN_mat12.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))

file.remove("ccw_presunuto.txt") ; file.remove("cw.txt")


## here, the files of the 1st and 2nd digitising persons are joined together

JN_YN_mat12<-rbind(JN_mat12,YN_mat12)


## Generalised Procrustes Analysis (GPA) of the entire dataset with sliding semilandmarks (using the minimum bending energy criterion)
# the file with information on the order of semilandmarks ("80_sliders.txt") is supplied by a user

pocetobjektu<-nrow(JN_YN_mat12) ; m<-2
pocetlandmarku<-ncol(JN_YN_mat12)/2
sliders<-read.table(file="80_sliders.txt")
h<-arrayspecs(JN_YN_mat12, pocetlandmarku, m, sep = NULL)
r<-gpagen(h, curves = sliders, PrinAxes = TRUE, ProcD = FALSE, Proj = TRUE)
e<-two.d.array(r$coords)
write.table(e, file = "aligned.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"))


## this part makes the segment configurations symmetric following Klingenberg et al. (2002, Evolution)

cm<-read.table(file="aligned.txt")

paroveLM<-c(1,1, 80,2, 79,3, 78,4, 77,5, 76,6, 75,7, 74,8, 73,9, 72,10, 71,11, 70,12, 69,13, 68,14, 67,15, 66,16, 65,17, 64,18, 63,19, 62,20, 61,21, 60,22, 59,23, 58,24, 57,25, 56,26, 55,27, 54,28, 53,29, 52,30, 51,31, 50,32, 49,33, 48,34, 47,35, 46,36, 45,37, 44,38, 43,39, 42,40, 41,41, 40,42, 39,43, 38,44, 37,45, 36,46, 35,47, 34,48, 33,49, 32,50, 31,51, 30,52, 29,53, 28,54, 27,55, 26,56, 25,57, 24,58, 23,59, 22,60, 21,61, 20,62, 19,63, 18,64, 17,65, 16,66, 15,67, 14,68, 13,69, 12,70, 11,71, 10,72, 9,73, 8,74, 7,75, 6,76, 5,77, 4,78, 3,79, 2,80) 

cm2<-cm
pocetobjektu<-nrow(cm) ; pocetsloupcu<-ncol(cm) ; pocetlandmarku<-pocetsloupcu/2
xLM<-length(paroveLM)
cmcm<-t(cm) ; cmcmcm<-as.data.frame(cmcm)

a<-1 ; b<-1
while (a<pocetobjektu+1) {
while (b<pocetsloupcu+1) {

cm2[a,b]<-(cm2[a,b]*-1) ; cm2[a,b+1]<-cm2[a,b+1]
b<-b+2 }
b<-1 ; a<-a+1 }

d<-1
while (d<(xLM+1)) {
p1<-paroveLM[d]
pmA<-cm2[,((p1*2)-1)] ; pmB<-cm2[,(p1*2)]
write(pmA, file = "pom.txt", ncolumns = pocetobjektu, append = TRUE)
write(pmB, file = "pom.txt", ncolumns = pocetobjektu, append = TRUE)
d<-d+2 }
dd<-read.table(file="pom.txt") ; dd2<-t(dd)
file.remove("pom.txt")

gg1<-append(cmcmcm,dd) ; gg4<-as.data.frame(gg1) ; gg5<-t(gg4) ; gg6<-as.data.frame(gg5)

write.table(gg6, file = "dvojnasobny-orr.txt", append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

mc<-read.table(file="dvojnasobny-orr.txt")
pocetobjektu2<-nrow(mc) ; m<-2

MC<-read.in("c:/R_prac/dvojnasobny-orr.txt",pocetlandmarku,m)

s<-c(pocetlandmarku,m,pocetobjektu2)
h<-array(data = MC , dim = s , dimnames = NULL)
r<-procGPA(h)
t<-r$rotated ; t2<-as.matrix(t)

m2<-1 ; i2<-nrow(t2) ; i2pol<-i2/2
while (m2<i2pol+1) {
kl1<-t2[m2] ; kl2<-t2[m2+i2pol]
kl3<-(kl1+kl2)/2
write(kl3, file = "pom2.txt", ncolumns = 1, append = TRUE)
m2<-m2+1 }

kl4<-read.table(file="pom2.txt")
kum<-matrix(data = NA, nrow = pocetobjektu, ncol = pocetsloupcu, byrow = FALSE)

kk<-1 ; ll<-1 ; mm<-1

while (mm<(nrow(kl4)-pocetlandmarku+1)) {
while (kk<pocetobjektu+1) {
while (ll<pocetsloupcu) {

gg<-kl4[mm,] ; gg2<-kl4[(mm+pocetlandmarku),] ; kum[kk,ll]<-gg ; kum[kk,ll+1]<-gg2
ll<-ll+2 ; mm<-mm+1 }
ll<-1 ; kk<-kk+1 ; mm<-mm+pocetlandmarku} }

file.remove("pom2.txt")
write.table(kum, file = "symetrizovany.txt", append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
file.remove("aligned.txt") ;  file.remove("dvojnasobny-orr.txt")


## this part will now conduct the Procrustes ANOVA decomposing the shape variation into multiple sources spanned by the independent factors

e<-read.table(file = "symetrizovany.txt")
pocetobjektu<-nrow(e)

faktory<-read.table(file="faktory.txt")

# region
reg<-faktory[,1] ; region<-as.vector(reg)
# locality nested in region
lok<-faktory[,2] ; lokalita<-as.vector(lok)
# plant nested in locality
ros<-faktory[,3] ; rostlina<-as.vector(ros)
# position as an independent fixed effect
poz<-faktory[,4] ; pozice<-as.vector(poz)
# position:region
pzr<-faktory[,7] ; pozice_region<-as.vector(pzr)
# position:region:locality
plk<-faktory[,8] ; pozice_region_lokalita<-as.vector(plk)
# position:region:locality:plant
plr<-faktory[,9] ; pozice_region_lokalita_rostlina<-as.vector(plr)
# error1 - digitising person
err_osoba<-read.table(file="err_osoba.txt")
osb<-err_osoba[,1] ; error_osoba<-as.vector(osb)
# error2 - digitising direction
err_digit<-read.table(file="err_digit.txt")
dig<-err_digit[,1] ; error_digit<-as.vector(dig)
# error1:error2
err_os_dig<-read.table(file="err_os_dig.txt")
eod<-err_os_dig[,1] ; error_os_dig<-as.vector(eod)


## this is the complete model; df, SS, MS and R2 are taken from here

h4<-procD.lm(e ~ pozice + region + lokalita + rostlina + pozice_region + pozice_region_lokalita +  pozice_region_lokalita_rostlina, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
summary(h4)

# secondary run for quantification of the measurement errors and the overall residuals

h4<-procD.lm(e ~ pozice_region_lokalita_rostlina + error_osoba + error_digit + error_os_dig, RRPP=T, int.first=T, seed = "random", iter=1)


## this part will create a number of random "region" factors by re-shuffling assignment of localities across regions 
# the model is then run with each randomly created "region" effect and the significance of the original effect may be assessed by comparison of the original SS with the distribution of random SS values

z<-1
while (z<pocetpermutaci+1) {

lok<-table(lokalita) ; lok2<-as.data.frame(lok)
x<-c("I1","I2","I3","M1","M2","M3","P1","P2","P3","V1","V2","V3")
X<-c("I","M","P","V")
ya<-sample(x,12,replace=F) ; # yb<-sample(x,12,replace=F)
# yc<-sample(x,12,replace=F) ; # yd<-sample(x,12,replace=F)

b<-1
while (b<length(x)+1) {
ya2<-ya[b] ; # yb2<-yb[b] ; # yc2<-yc[b] ; # yd2<-yd[b]
a<-1
while (a<(lok2[b,2]/4)+1) {
write(ya2, file = "rpa.txt", ncolumns = 1, append = TRUE)
a<-a+1 }
b<-b+1 }

fp_a<-read.table(file="rpa.txt")
fp_a2<-substr(t(fp_a),start=1,stop=1)
write(t(fp_a2), file = "region_perm_a.txt", ncolumns = 1, append = TRUE)
fp_a<-read.table(file="region_perm_a.txt")
# fp_b<-read.table(file="region_perm_b.txt")
# fp_c<-read.table(file="region_perm_c.txt")
# fp_d<-read.table(file="region_perm_d.txt")

rxp_a<-fp_a[,1] ; # rxp_b<-fp_b[,1] ; # rxp_c<-fp_c[,1] ; # rxp_d<-fp_d[,1]
rp_a<-as.vector(rxp_a) ; # rp_b<-as.vector(rxp_b)
# rp_c<-as.vector(rxp_c) ; # rp_d<-as.vector(rxp_d)
region_perm_FIN<-c(rp_a,rp_a,rp_a,rp_a)

# linear model with a random "region" factor

h5<-procD.lm(e ~ pozice + region_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=0)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals","Residuals"))
write.table(pok3$table$SS[2], file = "random_SSs_region.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

file.remove("region_perm_a.txt") ; file.remove("rpa.txt")
z<-z+1 }


## now, this part will create a number of random "locality" factors by re-shuffling assignment of plants across localities 

z<-1
while (z<pocetpermutaci+1) {

lok<-table(rostlina) ; lok2<-as.data.frame(lok)
x<-unique(rostlina)
X<-c("I1","I2","I3","M1","M2","M3","P1","P2","P3","V1","V2","V3")
ya<-sample(x,length(x),replace=F) ; # yb<-sample(x,length(x),replace=F)
# yc<-sample(x,length(x),replace=F) ; # yd<-sample(x,length(x),replace=F)

b<-1
while (b<length(x)+1) {
ya2<-ya[b] ; # yb2<-yb[b] ; # yc2<-yc[b] ; # yd2<-yd[b]
a<-1
while (a<(lok2[b,2]/4)+1) {
write(ya2, file = "rpa.txt", ncolumns = 1, append = TRUE)
a<-a+1 }
b<-b+1 }

fp_a<-read.table(file="rpa.txt")
fp_a2<-substr(t(fp_a),start=1,stop=2)
write(t(fp_a2), file = "lokalita_perm_a.txt", ncolumns = 1, append = TRUE)
fp_a3<-read.table(file="lokalita_perm_a.txt")
# fp_b<-read.table(file="lokalita_perm_b.txt")
# fp_c<-read.table(file="lokalita_perm_c.txt")
# fp_d<-read.table(file="lokalita_perm_d.txt")

fp_A2<-substr(t(fp_a),start=1,stop=1)
write(t(fp_A2), file = "region_perm_A.txt", ncolumns = 1, append = TRUE)
fp_A3<-read.table(file="region_perm_A.txt")

rxp_a<-fp_a3[,1] ; # rxp_b<-fp_b[,1] ; # rxp_c<-fp_c[,1] ; # rxp_d<-fp_d[,1]
rxp_A<-fp_A3[,1] ; # rxp_b<-fp_b[,1] ; # rxp_c<-fp_c[,1] ; # rxp_d<-fp_d[,1]

rp_a<-as.vector(rxp_a) ; # rp_b<-as.vector(rxp_b)
# rp_c<-as.vector(rxp_c) ; # rp_d<-as.vector(rxp_d)
rp_A<-as.vector(rxp_A)

lokalita_perm_FIN<-c(rp_a,rp_a,rp_a,rp_a)
region_perm_FIN<-c(rp_A,rp_A,rp_A,rp_A)

# linear model with a random "locality" factor

h5<-procD.lm(e ~ pozice + region + region_perm_FIN + lokalita_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals","Residuals","Residuals","Residuals"))
write.table(pok3$table$SS[4], file = "random_SSs_lokalita.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

file.remove("lokalita_perm_a.txt") ; file.remove("rpa.txt")
file.remove("region_perm_A.txt")
z<-z+1 }


## this part creates a number of random "plant" factors by re-shuffling assignment of segments across plants

z<-1
while (z<pocetpermutaci+1) {
x<-rostlina[1:(pocetobjektu/4)]
ya<-sample(x,length(x),replace=F)
y_a2<-substr(ya,start=1,stop=1) ; y_a3<-substr(ya,start=1,stop=2)

rostlina_perm_FIN<-c(ya,ya,ya,ya)
lokalita_perm_FIN<-c(y_a3,y_a3,y_a3,y_a3)
region_perm_FIN<-c(y_a2,y_a2,y_a2,y_a2)

# linear model with random "plant" factor

h5<-procD.lm(e ~ pozice + region + lokalita + region_perm_FIN + lokalita_perm_FIN + rostlina_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals","Residuals","Residuals","Residuals","Residuals","Residuals"))
write.table(pok3$table$SS[6], file = "random_SSs_rostlina.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

z<-z+1 }


## now, evaluation of "position"

z<-1
while (z<pocetpermutaci+1) {

x<-pozice[1:(pocetobjektu/4)]
ya<-sample(x,length(x),replace=F)
pozice_perm_FIN<-c(ya,ya,ya,ya)

# linear model with random "position" factor

h5<-procD.lm(e ~ pozice_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=0)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals"))
write.table(pok3$table$SS[1], file = "random_SSs_pozice.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
z<-z+1 }


## evaluation of "position:region"; randomization within "region"

z<-1
while (z<pocetpermutaci+1) {
x<-pozice_region[1:(pocetobjektu/4)]
y<-c("I","M","P","V")
a<-1
while (a<(length(y)+1)) {
b<-1 ; jj<-NULL
while (b<(length(x)+1)) {
if (region[b]==y[a]) { jj<-append(jj,x[b]) }
b<-b+1 }
jja<-sample(jj,length(jj),replace=F)
write(jja, file = "pozice_region_perm.txt", ncolumns = 1, append = TRUE)
a<-a+1 }

fp_a<-read.table(file="pozice_region_perm.txt") ; fp_a2<-t(fp_a)
fp_a3<-as.vector(fp_a2)
pozice_region_perm_FIN<-c(fp_a3,fp_a3,fp_a3,fp_a3)
pozice_perm_FIN<-substr(pozice_region_perm_FIN,start=1,stop=1)
region_perm_FIN<-substr(pozice_region_perm_FIN,start=3,stop=3)

# linear model with random "position:region" factor

h5<-procD.lm(e ~ pozice + region + lokalita + rostlina + pozice_perm_FIN + pozice_region_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals", "Residuals","Residuals","Residuals","Residuals","Residuals"))
write.table(pok3$table$SS[6], file = "random_SSs_pozice_region.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

file.remove("pozice_region_perm.txt")
z<-z+1 }


## evaluation of "position:region:locality"; randomization within "locality"

z<-1
while (z<pocetpermutaci+1) {
x<-pozice_region_lokalita[1:(pocetobjektu/4)]
y<-c("I1","I2","I3","M1","M2","M3","P1","P2","P3","V1","V2","V3")
a<-1
while (a<(length(y)+1)) {
b<-1 ; jj<-NULL
while (b<(length(x)+1)) {
if (lokalita[b]==y[a]) { jj<-append(jj,x[b]) }
b<-b+1 }
jja<-sample(jj,length(jj),replace=F)
write(jja, file = "pozice_region_lokalita_perm.txt", ncolumns = 1, append = TRUE)
a<-a+1 }

fp_a<-read.table(file="pozice_region_lokalita_perm.txt")
fp_a2<-t(fp_a) ; fp_a3<-as.vector(fp_a2)
pozice_region_lokalita_perm_FIN<-c(fp_a3,fp_a3,fp_a3,fp_a3)
pozice_perm_FIN<-substr(pozice_region_lokalita_perm_FIN,start=1,stop=1)
region_perm_FIN<-substr(pozice_region_lokalita_perm_FIN,start=3,stop=3)
lokalita_perm_FIN<-substr(pozice_region_lokalita_perm_FIN,start=3,stop=4)
pozice_region_perm_FIN<-substr(pozice_region_lokalita_perm_FIN,start=1,stop=3)

# linear model with a random "position:region:locality" factor

h5<-procD.lm(e ~ pozice + region + lokalita + rostlina + pozice_perm_FIN + pozice_region + pozice_region_perm_FIN + pozice_region_lokalita_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals", "Residuals","Residuals","Residuals","Residuals","Residuals","Residuals","Residuals"))
write.table(pok3$table$SS[8], file = "random_SSs_pozice_region_lokalita.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

file.remove("pozice_region_lokalita_perm.txt")
z<-z+1 }


## evaluation of "position:region:locality:plant"; randomization within "plant"

z<-1
while (z<pocetpermutaci+1) {
x<-pozice_region_lokalita_rostlina[1:(pocetobjektu/4)]
y<-unique(rostlina)
a<-1
while (a<(length(y)+1)) {
b<-1 ; jj<-NULL
while (b<(length(x)+1)) {
if (rostlina[b]==y[a]) { jj<-append(jj,x[b]) }
b<-b+1 }
jja<-sample(jj,length(jj),replace=F)
write(jja, file = "pozice_region_lokalita_rostlina_perm.txt", ncolumns = 1, append = TRUE)
a<-a+1 }

fp_a<-read.table(file="pozice_region_lokalita_rostlina_perm.txt")
fp_a2<-t(fp_a) ; fp_a3<-as.vector(fp_a2)
pozice_region_lokalita_rostlina_perm_FIN<-c(fp_a3,fp_a3,fp_a3,fp_a3)
pozice_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=1,stop=1)
region_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=3,stop=3)
lokalita_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=3,stop=4)
rostlina_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=3,stop=7)
pozice_region_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=1,stop=3)
pozice_region_lokalita_perm_FIN<-substr(pozice_region_lokalita_rostlina_perm_FIN,start=1,stop=4)

# linear model with random "position:region:locality:plant" factor

h5<-procD.lm(e ~ pozice + region + lokalita + rostlina + pozice_region + pozice_region_lokalita + pozice_perm_FIN + pozice_region_perm_FIN + pozice_region_lokalita_perm_FIN + pozice_region_lokalita_rostlina_perm_FIN, RRPP=F, SS.type = c("I"), int.first=F, seed = "random", iter=1)
pok3<-anova(h5, effect.type = c("F"), RRPP=F, error = c("Residuals", "Residuals","Residuals","Residuals","Residuals","Residuals","Residuals","Residuals","Residuals","Residuals"))
write.table(pok3$table$SS[10], file = "random_SSs_pozice_region_lokalita_rostlina.txt", append = TRUE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)

file.remove("pozice_region_lokalita_rostlina_perm.txt")
z<-z+1 }


