Ferramentas do usuário

Ferramentas do site


05_curso_antigo:alunos:trabalho_final:riguel

Riguel F Contente

riguel.jpg Doutorando em Oceanografia Biológica - IOUSP

Proposta principal

Mudanças alométricas relacionadas ao tamanho corpóreo podem afetar comparações de determinadas grandezas morfométricas entre populações. Lleonart et al. 2000 desenvolveram uma metodologia de remover esse efeito alométrico, a qual consiste na seguinte normalização: Z=Y(Tlo/Tl)^b, onde Y=valor da variável (e.g. tamanho da mandíbula), Tl= comprimento máximo do exemplar (em peixes, por exemplo, comprimento total), Tlo=representa um valor de referência (e.g. média ou mediana do comprimento total) através do qual todos as indivíduos são reduzidos. De um modo geral, o método remove do modelo o efeito do coeficiente 'a' ao qual se associa as informações de alometria. Pretendo realizar uma função que efetue essa normalização, produza gráficos de comparações de grandezas entre duas amostras (ex. boxplot, médias) e realize um teste de hipótese (teste-t ou anova) para avaliar possíveis diferenças significativas.

Referência

LLEONART J.; SALAT J.; TORRES G.J. Removing Allometric Effects of Body Size in Morphological Analysis Journal of Theoretical Biology, Volume 205, Number 1, 7 July 2000, pp. 85-93(9)

Comentários

Paulo
  1. Faça como o Fulano: formate a página no padrão que indicamos, po favor!
  2. Legal a proposta, e bem dimensionada. Apenas lembre que um teste t ou ANOVA tm várias premissas, e seria legal pelo indicar ao usuário se estas premissas parecem se cumprir com os dados.

Página de Ajuda

RemAllom                package:nenhum                R Documentation



Teste de alometria e comparação de grandezas morfométricas com remoção do efeito da alometria 

Description:

Realiza teste-t para avaliar possíveis desvios do crescimento isométrico, realiza ANOVA para testar se há relação entre variável resposta e preditora, remove potenciais efeitos de alometria através de uma normalização e, utilizando ANOVA, compara simultaneamente médias das estruturas morfométricas entre duas amostras/populações/espécies. Produz gráficos dos modelos alométricos e box-plots de comparação entre as amostras/populações/espécies. Exporta as matrizes com dados normalizados para cada população/espécie para uso em posteriores análises multivariadas.

Usage:
	
	RemAllom<-function(A,B,coef.A,coef.B,nomes,diag.allom=TRUE,TL0media=TRUE,diag.norm=TRUE,WT=TRUE)


Arguments:

	A: data.frame ou matriz da pop/espécie A. Estruturas morfométricas nas colunas e indivíduos nas linhas 
	B: data.frame ou matriz da pop/espécie B. Mesmo formato 
	coef.A: vetor numérico contendo os coeficientes b esperados para cada variável em A
	coef.B: vetor numérico contendo os coeficientes b esperados para cada variável em B
	nomes: objeto contendo a denominação para cada amostra/população/espécie
	diag.allom: lógico. Plotar gráficos diagnósticos para avaliação das premissas de regressão linear?
	TL0media: lógico. Reduzir os valores pela média? Se FALSE, será pela mediana
	diag.norm: lógico. Plotar gráficos exploratórios para avaliação de normalidade e homogeneidade das variâncias?
	WT: lógico. Exportar em .txt matrizes com valores normalizados?
	
Details:

Através do test-t, testa-se a hipótese nula de que o coeficiente b do modelo alométrico log-transformado (ln) é igual ao valor esperado de acordo com a dimensão da variável morfométrica em questão. Para comparação das variáveis entre duas populações/espécies, os dados são normalizados pela média ou mediana do comprimento do tamanho corpóreo, segundo a expressão:Z=Y(Tlo/Tl)^b, onde Y=valor da variável, Tl= comprimento corpóreo máximo do exemplar (ou outra medida de referência), Tlo= valor de referência, média ou mediana do comprimento total. O efeito da espécie/população sobre o tamanho médio de cada variável é posteriormente testado pela análise de variância.

Value:

Gráficos das regressões e box-plots de comparações serão gerados. Um data frame é gerado apresentando o resultado do teste-t para alometria e da análise de variância para avaliação de associação entre variável resposta e independente. Outro data frame exibe a tabela ANOVA de comparação entre pop/espécies. Dois arquivos .txt serão gerados no diretório de trabalho. 

Warning:

Atente se os dados exibem distribuição normal, e se as variâncias das variáveis em comparação são homogêneas. Efetue os testes diagnósticos para explorar os dados.
O designe da ANOVA é para amostras balanceadas, logo NAs (que serão automaticamente removidos) podem tornar a análise inválida a subsequentes interpretações. 
Note que não há dimensão nos box-plot, pois Z é adimensional

Author(s):

Riguel Feltrin Contente

riguel@io.usp.br


References:

Slope of a regression line t-test - http://en.wikipedia.org/wiki/T-test
 
LLEONART J.; SALAT J.; TORRES G.J. Removing Allometric Effects of Body Size in Morphological Analysis. Journal of Theoretical Biology, Volume 205, Number 1, 7 July 2000, pp. 

85-93(9)

Examples:

data(iris) ## Abrir "iris"
setosa1=iris[iris$Species=="setosa",] ## Selecionar apenas os vetores com as variáveis morfométricas da sépala e pétala da espécie 'setosa' 
setosa<-setosa[,-5]
head(setosa)
virginica1=iris[iris$Species=="virginica",] ## Selecionar apenas os vetores com as variáveis morfométricas da sépala e pétala da espécie 'virginica'
virginica<-virginica[,-5]
head(virginica)
setosa.b<-c(1,1,1,1) ## Criar um vetor numérico com os coeficientes esperados de acordo com a dimensão da variável 
virginica.b<-c(1,1,1,1)
nomes<-c("Setosa","Virginica") ## Criar objeto com o nome das espécies
RemAllom(setosa,virginica,setosa.b,virginica.b,nomes,diag.allom=F,TL0media=T,diag.norm=F,WT=T) ## Rodar a função na qual a normalização será efetuada pela média da sépala

Código da Função


RemAllom<-function(A,B,coef.A,coef.B,nomes,diag.allom=T,TL0media=T,diag.norm=T,WT=T)
{
##Matriz A
##Remoção das linhas com NA, transformação ln e calculo da media e mediana
	x.1<-na.omit(A)
	x<-log(x.1)
	MM<- matrix(rep(NA,dim(x)[2]*2),nrow=2)
	for(j in 1:dim(x)[2])
	{
		MM[1,j]=mean(x[,j])
		MM[2,j]=median(x[,j])
	}
##Remoção das linhas com NA e calculo da media e mediana
	x.3<-na.omit(A)
	MM2<- matrix(rep(NA,dim(x.3)[2]*2),nrow=2)
	for(j in 1:dim(x.3)[2])
	{
		MM2[1,j]=mean(x.3[,j])
		MM2[2,j]=median(x.3[,j])
	}
	mm=MM2
	med=MM
##Teste de hipótese (test-t) para crescimento alométrico 
	Stat<- matrix(rep(NA,(dim(x)[2]*8)-8),nrow=8)
	for(j in 2:dim(x)[2])
	{
		som.q<-sum((x[,1]-med[1,j])^2)
		mod.l<-lm(x[,j]~x[,1],data=x)
		gl<-length(x[,1])-2
		test.coef<-((coef(mod.l)[2]-coef.A[j])*sqrt(gl))/(sqrt(sum(residuals(mod.l)^2)/som.q))
		p.1<-2*pt(test.coef,df=gl,lower.tail=F)
		Stat[1,c(j-1)]=round(test.coef,3)
		Stat[2,c(j-1)]=round(p.1,5)
		Stat[3,c(j-1)]=round(coef(mod.l)[2],3)
		Stat[4,c(j-1)]=coef.A[j]
		Stat[6,c(j-1)]<-round(summary(mod.l)$r.squared,3)
		Stat[7,c(j-1)]<-round(anova(mod.l)$"F value"[1],3)
		Stat[8,c(j-1)]<-round(anova(mod.l)$"Pr(>F)"[1],7)
		StatA<-data.frame(Stat)
	}
	for(i in 1:dim(StatA)[2])
	{
		if(StatA[2,i]>=0.05)
			{StatA[5,i]=c("isom")}
		else{
			if(StatA[2,i]<0.05&StatA[3,i]>StatA[4,i])
				{StatA[5,i]=c("alom+")}
			else{
				if(StatA[2,i]<0.05&StatA[3,i]<StatA[4,i])
					{StatA[5,i]=c("alom-")
		}}}
	}
	for(k in 1:dim(StatA)[2])
	{
		if(StatA[2,k]>=0.05)
			{StatA[2,k]=c(">0.05")}
		else{
			if(StatA[2,k]<0.05&StatA[2,k]>0.0001)
				{StatA[2,k]=StatA[2,k]}	
			else{		
				if(StatA[2,k]<0.0001)
					{StatA[2,k]=c("<0.0001")
		}}}
 		if(StatA[8,k]>=0.05)
			{StatA[8,k]=c(">0.05")}
		else{
			if(StatA[8,k]<0.05&StatA[8,k]>0.0001)
				{StatA[8,k]=StatA[8,k]}	
			else{		
				if(StatA[8,k]<0.0001)
					{StatA[8,k]=c("<0.0001")
		}}}
	}
	StatB<- matrix(rep(NA,(dim(x)[2]*5)-5),nrow=5)
	for(j in 2:dim(x)[2])
	{
		som.q<-sum((x[,1]-med[1,j])^2)
		mod.l<-lm(x[,j]~x[,1],data=x)
		gl<-length(x[,1])-2
		test.coef<-((coef(mod.l)[2]-coef.A[j])*sqrt(gl))/(sqrt(sum(residuals(mod.l)^2)/som.q))
		p.1<-2*pt(test.coef,df=gl,lower.tail=F)
		StatB[1,c(j-1)]=round(test.coef,3)
		StatB[2,c(j-1)]=round(p.1,5)
		StatB[3,c(j-1)]=coef(mod.l)[2]
		StatB[4,c(j-1)]=coef.A[j]
	}
##Gráficos diagnósticos para o Modelo Gaussiano 
	if(diag.allom==T)
	{
		for(j in 2:dim(x)[2])
		{
			win.graph(width = 9, height = 9)
			par(mfrow=c(2,2))
			mlw<-lm(x[,j]~x[,1],data=x)
			plot(mlw,main=paste("ln",names(x[j]),nomes[1]))
		}
		par(mfrow=c(1,1))
	}
##Regressões 
	g=1
	for(i in 1:dim(x)[2])
	{
		g=g+1
		if(g<=dim(x)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			mlw<-lm(x[,g]~x[,1],data=x)
			gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
			mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
			abline(mlw,col=2)
			g=g+1
			if(g<=dim(x)[2]){
				mlw<-lm(x[,g]~x[,1],data=x)
				gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
				mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
				abline(mlw,col=2)
				g = g+1
				if(g<=dim(x)[2]){
					mlw<-lm(x[,g]~x[,1],data=x)
					gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
					mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
					abline(mlw,col=2)
					g = g+1
					if(g<=dim(x)[2]){
						mlw<-lm(x[,g]~x[,1],data=x)
						gra<-plot(x[,g]~x[,1],data=x,main=nomes[1],ylab=paste("ln",names(x[g])),xlab=paste("ln",names(x[1])))
						mtext(StatA[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
						abline(mlw,col=2)
		}}}}
	}
##Dados normalizados pela média ou mediana
	R.all<-matrix(rep(NA,dim(x.3)[1]*dim(x.3)[2]),nrow=dim(x.3)[1],ncol=dim(x.3)[2])
	for(l in 1:dim(x.3)[1])		
	{
		for(j in 2:dim(x.3)[2])
		{
			if(TL0media==T)
				{raz<-x.3[l,j]*((mm[1,1]/x.3[l,1])^StatB[3,j-1])
				R.all[l,j-1]<-raz}
			else{
				{raz<-x.3[l,j]*((mm[2,1]/x.3[l,1])^StatB[3,j-1])
				R.all[l,j-1]<-raz}}
		}
	}
##Resultados - dados matriz A	
	nomes.0<-c(nomes[1],"","","","","","","")
	nomes.1<-c("-->TESTE-T","p/ avaliar","alometria","","","-->ANOVA","p/ avaliar","sig. regressão")
	nomes.2<-c("t","Pr(>|t|)","b obs","b exp","cresc","r2","F","Pr(>F)")
	StatC<-data.frame(nomes.0,nomes.1,nomes.2,StatA)
	nomes.componentes<-c("sp/pop.","testes","comp.")
	nomes.col<-matrix(rep(NA,(dim(x)[2])-1),nrow=1)
	for(j in 2:dim(x)[2])
	{
		nomes.col[1,j-1]<-names(x[j])
	}
	MF<-R.all[,-c(dim(R.all))]
	colnames(MF)<-nomes.col
	names(StatC)<-c(nomes.componentes,nomes.col)

##Teste de normalidade- dados matriz A
	if(diag.norm==T)
		{x11()
		par(tcl=0.3)
		par(mfrow=c(2,2),pty="s")
		devAskNewPage(TRUE)
		for(m in 1:dim(MF)[2])
		{
			qqnorm(MF[,m],main=paste(nomes[1],nomes.col[1,m]))
			qqline(MF[,m],col=2)
			hist(MF[,m],prob=T,main="",xlab=nomes.col[1,m])
			par(new=T)
			curve(expr=dnorm(x,mean(MF[,m]),sd(MF[,m])),main="",col=2,ann=F,axes=F)
		}
		par(mfrow=c(1,1))
	}
	devAskNewPage(F)

##Matriz B
##Remoção das linhas com NA, transformação ln e calculo da media e mediana
	y.1<-na.omit(B)
	y<-log(y.1)
	MMy<- matrix(rep(NA,dim(y)[2]*2),nrow=2)
	for(j in 1:dim(y)[2])
	{
		MMy[1,j]=mean(y[,j])
		MMy[2,j]=median(y[,j])
	}
##Remoção das linhas com NA e calculo da media e mediana
	y.3<-na.omit(B)
	MM2y<- matrix(rep(NA,dim(y.3)[2]*2),nrow=2)
	for(j in 1:dim(y.3)[2])
	{
		MM2y[1,j]=mean(y.3[,j])
		MM2y[2,j]=median(y.3[,j])
   	}
##Teste de hipótese (test-t) para crescimento alométrico 
	mmy=MM2y
	medy=MMy
	Staty<- matrix(rep(NA,(dim(y)[2]*8)-8),nrow=8)
	for(j in 2:dim(y)[2])
	{
		som.qy<-sum((y[,1]-medy[1,j])^2)
		mod.ly<-lm(y[,j]~y[,1],data=y)
		gly<-length(y[,1])-2
		test.coefy<-((coef(mod.ly)[2]-coef.B[j])*sqrt(gly))/(sqrt(sum(residuals(mod.ly)^2)/som.qy))
		p.1y<-2*pt(test.coefy,df=gly,lower.tail=F)
		Staty[1,c(j-1)]=round(test.coefy,3)
		Staty[2,c(j-1)]=round(p.1y,5)
		Staty[3,c(j-1)]=round(coef(mod.ly)[2],3)
		Staty[4,c(j-1)]=coef.B[j]
		Staty[6,c(j-1)]<-round(summary(mod.ly)$r.squared,3)
		Staty[7,c(j-1)]<-round(anova(mod.ly)$"F value"[1],3)
		Staty[8,c(j-1)]<-round(anova(mod.ly)$"Pr(>F)"[1],7)
		StatAy<-data.frame(Staty)
	}
	for(i in 1:dim(StatAy)[2])
	{
		if(StatAy[2,i]>=0.05)
			{StatAy[5,i]=c("isom")}
		else{
			if(StatAy[2,i]<0.05&StatAy[3,i]>StatAy[4,i])
				{StatAy[5,i]=c("alom+")}
			else{
				if(StatAy[2,i]<0.05&StatAy[3,i]<StatAy[4,i])
					{StatAy[5,i]=c("alom-")
		}}}
	}
	for(k in 1:dim(StatAy)[2])
	{
		if(StatAy[2,k]>=0.05)
			{StatAy[2,k]=c(">0.05")}
		else{
			if(StatAy[2,k]<0.0499999&StatA[2,k]>0.0001)
				{StatAy[2,k]=StatAy[2,k]}	
			else{		
				if(StatAy[2,k]<0.0001)
					{StatAy[2,k]=c("<0.0001")
		}}}
		if(StatAy[8,k]>=0.05)
			{StatAy[8,k]=c(">0.05")}
		else{
			if(StatAy[8,k]<0.0499999&StatA[8,k]>0.0001)
				{StatAy[8,k]=StatAy[8,k]}	
			else{		
				if(StatAy[8,k]<0.0001)
					{StatAy[8,k]=c("<0.0001")
		}}}
	}
	StatBy<- matrix(rep(NA,(dim(y)[2]*5)-5),nrow=5)
	for(j in 2:dim(y)[2])
	{
		som.qy<-sum((y[,1]-medy[1,j])^2)
		mod.ly<-lm(y[,j]~y[,1],data=y)
		gly<-length(y[,1])-2
		test.coefy<-((coef(mod.ly)[2]-coef.B[j])*sqrt(gly))/(sqrt(sum(residuals(mod.ly)^2)/som.qy))
		p.1y<-2*pt(test.coefy,df=gly,lower.tail=F)
		StatBy[1,c(j-1)]=round(test.coefy,3)
		StatBy[2,c(j-1)]=round(p.1y,5)
		StatBy[3,c(j-1)]=coef(mod.ly)[2]
		StatBy[4,c(j-1)]=coef.B[j]
	}
##Gráficos diagnósticos para o Modelo Gaussiano 
	if(diag.allom==T)
	{
		for(j in 2:dim(y)[2])
		{
			win.graph(width=9,height=9)
			par(mfrow=c(2,2))
			mlwy<-lm(y[,j]~y[,1],data=y)
			plot(mlwy,main=paste("ln",names(x[j]),nomes[2]))
		}
		par(mfrow=c(1,1))
	}
##Regressões		
	g=1
	for(i in 1:dim(y)[2])
	{
		g=g+1
		if(g<=dim(y)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			mlwy<-lm(y[,g]~y[,1],data=y)
			grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
			mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
			abline(mlwy,col=2)
			g=g+1
			if(g<=dim(y)[2]){
				mlwy<-lm(y[,g]~y[,1],data=y)
				grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
				mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
				abline(mlwy,col=2)
				g = g+1
				if(g<=dim(y)[2]){
					mlwy<-lm(y[,g]~y[,1],data=y)
					grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
					mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
					abline(mlwy,col=2)
					g = g+1
					if(g<=dim(y)[2]){
						mlwy<-lm(y[,g]~y[,1],data=y)
						grayy<-plot(y[,g]~y[,1],data=y,main=nomes[2],ylab=paste("ln",names(y[g])),xlab=paste("ln",names(y[1])))
						mtext(StatAy[5,g-1],side=3,cex=0.8,line=0.75,las=0,font=2,adj=0)
						abline(mlwy,col=2)
		}}}}
	}
##Dados normalizados pela média ou mediana
	R.ally<-matrix(rep(NA,dim(y.3)[1]*dim(y.3)[2]),nrow=dim(y.3)[1],ncol=dim(y.3)[2])
	for(l in 1:dim(y.3)[1])		
	{
		for(j in 2:dim(y.3)[2])
		{
			if(TL0media==T)
				{razy<-y.3[l,j]*((mmy[1,1]/y.3[l,1])^StatBy[3,j-1])
				R.ally[l,j-1]<-razy}
			else{
				{razy<-y.3[l,j]*((mmy[2,1]/y.3[l,1])^StatBy[3,j-1])
				R.ally[l,j-1]<-razy}}
		}
	}	
	nomes.0y<-c(nomes[2],"","","","","","","")
	nomes.1y<-c("-->TESTE-T","p/ avaliar","alometria","","","-->ANOVA","p/ avaliar","sig. regressão")
	nomes.2y<-c("t","Pr(>|t|)","b obs","b exp","cresc","r2","F","Pr(>F)")
	StatCy<-data.frame(nomes.0y,nomes.1y,nomes.2y,StatAy)
	nomes.componentesy<-c("sp/pop.","testes","comp.")
	nomes.coly<-matrix(rep(NA,(dim(y)[2])-1),nrow=1)
	for(j in 2:dim(y)[2])
	{
		nomes.coly[1,j-1]<-names(y[j])
	}
	MFy<-R.ally[,-c(dim(R.ally))]
	colnames(MFy)<-nomes.coly
	names(StatCy)<-c(nomes.componentesy,nomes.coly)

##Teste de normalidade- dados matriz B
	if(diag.norm==T)
		{x11()
		par(tcl=0.3)
		par(mfrow=c(2,2),pty="s")
		devAskNewPage(TRUE)
		for(m in 1:dim(MFy)[2])
		{
			qqnorm(MFy[,m],main=paste(nomes[2],nomes.coly[1,m]))
			qqline(MFy[,m],col=2)
			hist(MFy[,m],prob=T,main="",xlab=nomes.coly[1,m])
			par(new=T)
			curve(expr=dnorm(x,mean(MFy[,m]),sd(MFy[,m])),main="",col=2,ann=F,axes=F)
		}
		par(mfrow=c(1,1))}
		devAskNewPage(F)

##Geração dos dados normalizados e análises comparativas
##Geração do fator
	MFA=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			MFA[i,j]=nomes[1]
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			MFA[i+dim(MF)[1],j]=nomes[2]
			as.factor(MFA)
			MFB<-MFA[,1]
	}}
##Matriz p/ as análises
	mtx=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mtx[i,j]=MF[i,j]
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mtx[i+dim(MF)[1],j]=MFy[i,j]
	}}	
##Exportação das matrizes com dados normalizados
	
	if(WT==T)
		{write.table(MF,file="MFA.txt",row.names=T,col.names=T,sep=",")
		write.table(MFy,file="MFB.txt",row.names=T,col.names=T,sep=",")
	}
	
##Box-plots
	g=0
	for(i in 1:dim(mtx)[2])
	{
		g=g+1
		if(g<=dim(mtx)[2]){
			x11()
			par(mfrow=c(2,2))
			par(tcl=0.3)
			par(mgp=c(1.8,0.5,0))
			boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
			g=g+1
			if(g<=dim(mtx)[2]){
				boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
				g = g+1
				if(g<=dim(mtx)[2]){
					boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
					g = g+1
					if(g<=dim(mtx)[2]){
						boxplot(mtx[,g]~MFB,main=nomes.col[g],ylab=names(mtx[,g]))
		}}}}
	}

##ANOVA
	
##SS Total
	media.geral=matrix(rep(NA,dim(MF)[2]),nrow=1)
	for(i in 1:dim(MF)[2])
	{
		media.geral[1,i]=mean(c(MF[,i],MFy[,i]))
	}
	mtx1=matrix(rep(NA,(dim(MF)[2])*(dim(MF)[1]*2)),nrow=dim(MF)[1]*2,ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mtx1[i,j]=(MF[i,j]-media.geral[1,j])^2
	}}
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mtx1[i+dim(MF)[1],j]=(MFy[i,j]-media.geral[1,j])^2
	}}	
	ss.total=(apply(mtx1,2,sum))

##Somatória Quadráticas Intra grupo
	mediaA.q<-matrix(rep(NA,(dim(MF)[1]*dim(MF)[2])),nrow=dim(MF)[1],ncol=dim(MF)[2])
	for(i in 1:dim(MF)[1])
	{
		for(j in 1:dim(MF)[2])
		{	
			mediaA.q[i,j]<-(MF[i,j]-mean(c(MF[,j])))^2
	}}
	ss.mediaA.q=apply(mediaA.q,2,sum)
	mediaB.q<-matrix(rep(NA,(dim(MFy)[1]*dim(MFy)[2])),nrow=dim(MFy)[1],ncol=dim(MFy)[2])
	for(i in 1:dim(MFy)[1])
	{
		for(j in 1:dim(MFy)[2])
		{	
			mediaB.q[i,j]<-(MFy[i,j]-mean(c(MFy[,j])))^2
	}}
	ss.mediaB.q=apply(mediaB.q,2,sum)
	ss.intra=ss.mediaA.q+ss.mediaB.q
##Somatória Quadráticas Entre grupos
	mediaA<-apply(MF,2,mean)
	mediaB<-apply(MFy,2,mean)
	ss.entre=dim(MF)[1]*((mediaA-media.geral)^2+(mediaB-media.geral)^2)
##Cálculo do F
	GL1=1
	GL2=(dim(mtx)[1]-2)
	ms.entre=ss.entre/GL1
	ms.intra=ss.intra/GL2
	F.anova=ms.entre/ms.intra
	
	format.pval(c(0.1, 0.001, 1e-27))
	
##Cálculo do P
	pf.anova=matrix(rep(NA,dim(MF)[2]*1),ncol=dim(MF)[2],nrow=1)
	for(i in 1:dim(MF)[2])
	{		
		pf.anova[1,i]=pf(F.anova[1,i],GL1,GL2,lower.tail=FALSE)
	}
##Tabela ANOVA
	PPato=round(c(pf.anova),3)
	SSE=round(c(ms.entre),3)
	MSI=round(c(ms.intra),3)
	FI=round(c(F.anova),3)
	SSI=round(c(ss.intra),3)
	GL1m=rep(1,dim(MF)[2])
	GL2m=rep(GL2,dim(MF)[2])
	GL1.r=round(GL1m,1)
	GL2.r=round(GL2m,1)	
	Tabela.ANOVA<-matrix(c(SSE,SSI,GL1.r,GL2.r,SSE,MSI,FI,PPato),ncol=dim(MF)[2],nrow=8,byrow=T)
	colnames(Tabela.ANOVA)=nomes.col
	rownames(Tabela.ANOVA)=c("Sq.entre","Sq.res","gl1","gl2","Mq.entre","Mq.res","F","Pr(>F)")
##Resultado 
	cat("\t\nResultado teste-t - avaliação de desvios do crescimento isométrico \nResultado ANOVA - teste de associação entre as variáveis resposta e independente \nTabela ANOVA\n")
	return(StatC,StatCy,Tabela.ANOVA)
	
}	

Arquivo da Função

05_curso_antigo/alunos/trabalho_final/riguel.txt · Última modificação: 2020/08/12 06:04 (edição externa)