anova fatorial e problemas em achar funções

Boa tarde porfessor Marcos,
Estou escrevendo um artigo e preciso comparar riqueza e abundância de macroinvertebrados bentônicos que coltei em 3 rios com os dados das amostras coletados nesses três rios em outros dois períodos diferentes. Ao meu ver eu conseguiria fazer essa análise com uma anova fatorial, mas esbarrei com duas duvidas:
A primeira é que eu fiz o teste de Levine e ele apontou homogeneidade da riqueza em relação aos rios, mas esse pressuposto não foi cumprido em relação aos períodos, então pensei em fazer um análise GLS, até rodei ela mas o gráfico ficou muito diferente, nada parecido com os que vimos nas aulas. Devo fazer desse modo mesmo ou devo fazer diferente?
A segunda é que na hora de gerar os gráficos, o meu RStudio não reconhece as funções lineplot.CI e sciplot, como eu faço para corrigir isso?

Montei minha tabela dessa forma:

corrego periodo riqueza abundancia
1m pre 163 11052
bale pre 154 5252
nsp pre 141 2515
1m pos1 52 7007
bale pos1 95 6709
nsp pos1 43 1064
1m pos2 32 200
bale pos2 28 198
nsp pos2 35 790

Meu script ficou dessa forma. mas eu acho que tem alguns erros ai, será que você poderia me ajudar?

###########################################
# Resumo Simpósio Bacia Hidrográfica do rio São Francisco #
# Teste ANOVA #
# Comparando riqueza entre três córregos em diferentes períodos #
# Moana Rothe - Neves #
###########################################

#Definindo o diretório:
setwd("D:/Dados do Pc/Documents/III SBHSF/análises R")
dir()

#Carregando os dados:
dados<-read.table("riqueza e abundancia.txt", h=T)
dados
summary(dados)

####################

#Para se fazer uma ANOVA é preciso testar alguns pressupostos:
#Homogeneidade das variâncias
#Normalidade dos resíduos (esse é feito depois da ANOVA)

####################

#Testando os pressupostos da ANOVA:

#Homogeneidade:
library(car)
leveneTest(dados$riqueza ~ dados$corrego)
leveneTest(dados$riqueza ~ dados$periodo)

#Houve quebra do pressuposto da homogeneidade da riqueza em relação ao período
#Sendo assim, farei uma análise de GLS (Generalized Least Square)
#Ela considera a variância dos residuos separadamente.

#Para isso vou carregar o pacote nlme
library(nlme)

test1<-gls(riqueza~periodo, weights=varIdent(form=~3|periodo), data=dados)
plot(test1)
anova(test1)
summary(test1)

#Agora vou vizualizar através de um gráfico:
#Lembrando que por deffaut a barra de variação é o erro padrão, mas é possivel mudar

library(sciplot)

tiff('Fig.1.tif', w = 2400, h = 1800, compression = 'lzw', res = 300)
par(bty = "l", mar=c(4, 4, 2, 2), mgp = c(2.5, 1, 0))
lineplot.CI(factor(periodo, levels=c('pre', 'pos1', 'pos2')), riqueza, data=dados, xlab = "Córregos", ylab = "Riqueza", type = 'p')
dev.off()

#Criando o modelo de ANOVA fatorial:
resultado.anova<-aov(dados$riqueza ~ dados$corrego * dados$periodo)

#Normalidade dos resíduos:
qqnorm(resultado.anova$residuals)
qqline(resultado.anova$residuals, lty=2)

#Vamos avaliar o resultado:
summary(resultado.anova)

#Vamos entender o resultado:

lineplot.CI(dados$riqueza, dados$corrego, dados$periodo, las=1, xlab="corrego", ylab="riqueza")

lineplot.CI(dados$Sexo, dados$atividade, dados$Especie, las=1, xlab="Sexo", ylab="atividade")