#################################### ## COMEÇAR POR CORRER ESTA FUNÇÃO ## #################################### color_biplot<-function(data, axes=1:2, cex.ind=1, cex.var=1, labels.ind=1:nrow(X), col.ind=rep("black",nrow(data)), col.var=rep("red",ncol(data))){ data<-scale(data,scale=FALSE) M<-nrow(data) q<-ncol(data) X<-data N <-nrow(X) p <-ncol(X) X.svd <- svd(X) X.pca<-prcomp(X) U <- X.svd$u D <- diag(X.svd$d) V<- X.svd$v scores <- U %*% D sdev<-diag(D)/sqrt(N-1) ### PC's standard deviations = X.pca$sdev V.aux<-V G<-data %*% V.aux %*% solve(D)*sqrt(N-1) H=V.aux %*% D /max(abs(V.aux %*% D)[,axes])*max(abs(G)[,axes])/2 plot(G[,axes],cex=.001,asp=TRUE,xlab="PC1",ylab="PC2") abline(h=0,lty=2,col="gray") abline(v=0,lty=2,col="gray") text(G[,axes[1]],G[,axes[2]],col=col.ind,cex=cex.ind,labels.ind) arrows(rep(0,p),rep(0,p),H[,axes[1]],H[,axes[2]],col=col.var,length=.1) text(H[,axes[1]]*1.3,H[,axes[2]]*1.3,col=col.var,cex=cex.var,labels=colnames(X)) } #################################### ### 1 aula (15 de junho de 2026) ### #################################### # Resolução com comentários adicionais do exercício 7 ## Para lermos os dados sobre relativos às 9 culturas em 20 concelhos ## do distrito de santarem contidos no ficheiro dadosMulti.rdata ## (exercícios do Prof. Cadima), podemos executar a instrução abaixo load("/Users/pedrosilva/Documents/ISA25_26/MMA/Aulas/aulas/dadosMulti.rdata") ## Obs: não é necessário indicar o caminho das diretorias (apenas o ## nome do ficheiro) se indicarmos ao RStudio a diretoria que contém ## o ficheiro como diretoria de trabalho ls() ### mostra os objetos que foram carregados para o ambiente do RStudio X=as.matrix(santarem) ## converte o dataset numa data matrix View(X) ## para mostrar apenas as primeiras 6 linhas de X head(X) ## para mostrar apenas as últimas 6 linhas de X tail(X) ## a) ## Construção da ACP sobre os dados originais, ## isto é, sobre a matriz de covariâncias: X_ACP<-prcomp(X) ## é o mesmo que fazer X_ACP<-prcomp(X,scale=FALSE) X_ACP ### devolve um objeto com várias componentes: # desvios-padrão das PCs = sqrt(variancias) = sqrt(lambdas) sdev<-X_ACP$sdev sdev ## matriz dos loadings, isto é, ## dos coeficientes de cada uma ## das PCs como CL das variavéis originais loadings<-X_ACP$rotation loadings # matriz dos scores, isto é, ## das coordenadas dos 20 indivíduos (concelhos) # nas novas variáveis sintéticas, PC1,...,PC9 scores<-X_ACP$x scores ## vetor das médias das 9 culturas, isto é, ## baricentro ou centro de gravidade da núvem de 20 pts em |R^9 pois ## cada um dos 20 concelhos está associado a um vetor com 9 componentes ## que representa as produções, em t/ha, das 9 culturas X_ACP$center colMeans(X) ### idem ## EXTRAS ### as variancias de PC1,..., PC9 (lambda1,...,lambda9) podem ### ser obtidas de 3 formas distintas: sdev^2 diag(var(scores)) S <- 1/(nrow(X)-1)*t(Xc) %*% Xc ## matriz de covariâncias de X - slide 19 eigen(S)$values ### valores próprios da matriz de covariâncias de X ## construção alternativa da matriz dos scores (ver o slide 34) Xc <- scale(X,scale=F) ## matriz centrada dos dados scores_bis <- Xc %*% loadings scores_bis ## deve ser igual à matriz dos scores acima ## i) summary(X_ACP) # variancias explicadas por cada PC e variancias acumuladas ### PC1 explica 94.2% (aprox) da variancia total ### a PC1+PC2 98.4% (aprox), etc... ### a qualidade da redução para duas dimensões é excelente ## A matriz de covariâncias mostra que o vetor de observações ## associado à cultura da batata var(X[,"batata"]) ## apresenta uma variância muito superior às restantes culturas diag(cov(X)) ## correspondendo, alías, a mais de 93% da variancia total var(X[,9])/sum(diag(cov(X))) ## Logo a PC1 que é uma combinação linear (CL) das 9 variáveis ## que maximiza a variância explicada não pode ## explicar menos que a CL definida só pela cultura batata, isto é, ## 0 * trigo + 0 * milho + ... + 0 * grao + 1 * batata ## que já explica 93% ## O facto da cultura da batata ter uma variância com um ## peso muito superior às restantes culturas atribui uma ## importância excessiva a esta variável enviasando a PCA, ## pelo que a PCA sobre os dados originais não é adequada ## (apesar das variáveis estarem todas expressas nas mesmas unidades). ## ii) scores[,1:2] ## coordenadas dos 20 concelhos relativos a PC1 e PC2 plot(scores[,1:2],asp=TRUE,col="red",pch=16) ## asp=TRUE -> keep the aspect ratio abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy ind=identify(scores[,1:2]) ## fica a aguardar que o utilizador aponte com o ponteiro para ## os pontos que prete nde identificar (para terminar premir ESC) rownames(X) ## nomes das linhas do data matrix X (concelhos - indiv) colnames(X) ## nomes das colunas do data matrix X (culturas - variav) rownames(X)[ind] # nomes das linhas do data matrix X que foram identificados ## escreve o nome por baixo (pos=1) do corresponde ponto no scatterplot text(scores[ind,1],scores[ind,2],labels=rownames(X)[ind],pos=1,cex=.5) ################################################################ ### 2 aula (16 de junho de 2026) - continuação do exercício 7 ## ################################################################ ## iii) ## Correlações entre as variáveis x1(trigo),..,x9 (batata) e as PCs ## (arredondadas às 5 casas décimais) Correl<-round(cor(X,scores),5) Correl ## Cálculo alternativo das correlações anteriores usando a ## fórmula do slide 37 ## cor(x_j,PCk)=sqrt(lambda_k) * v_jk / s_j ## x9=X[,9] -> batata ## x2=X[,2] -> milho ## v_jk: loading da variável x_j na PCk ## elemento da linha j e coluna k da matriz dos loadings ## Então a cor(x9,CP1) vem dada por: sdev[1]*loadings[9,1]/sd(X[,9]) ### ou sdev[1]*loadings[9,1]/sd(X[,"batata"]) # comparar com Correl[9,1] ## cor(x2,CP1) - milho X_ACP$sdev[1]*loadings[2,1]/sd(X[,"milho"]) Correl[2,1] ## cor(x6,CP1) - fava X_ACP$sdev[1]*loadings[6,1]/sd(X[,"fava"]) Correl[6,1] ## iv ## dada a correlação quase perfeita entre variável PC1 ## e a variável batata podemos ## interpretar o primeiro eixo como estando associado à produção ## da batata e nomeá-lo como tal. ### Analogamente o 2º eixo pode ser considerado como o eixo relativo ### à produção de milho dada a elevada correlação as variáveis ### milho e PC2 (0.93958) ## os loadings confirmam que a PC1 é CL quase exclusivamente ## obtida a partir da variável batata com um coeficiente de 0.9966 ## e coeficientes quase nulos relativamente às restantes variáveis. round(loadings[,1],4) ### idem relativamente ao milho round(loadings[,2],4) ## A análise do biplot ajudará a consolidar estas conclusões... ## v biplot(X_ACP) abline(v=0,lty=2) abline(h=0,lty=2) ## Por omissão o comando biplot efetua o biplot focado nas correlações ## e covariâncias (slide 61 e seguintes), sendo a qualidade da redução ## quase perfeita, uma vez que a variância explicada pelas componentes ## PC1 e PC2 é superior a 98%. ## Podemos assim aplicar as seguintes interpretações do slide 64: ## -------------------------------------------------------------- ### a correlação entre a variável batata e a PC1 é dada pelo cosseno do ### ângulo formado pelo marcador da variável batata e o eixo ### horizontal sendo o seu valor muito próximo de 1, em coerência ### com o valor dessa correlação calculado anteriormente. ### a correlação entre as variáveis batata e milho é quase nula, ### uma vez os respetivos marcadores são quase ortogonais entre si, ### isto é, o cosseno o ângulo é quase nulo. ### os comprimentos dos marcadores das variáveis são proporcionais ### aos desvios-padrão dessas variáveis constatando-se novamente ### que a variável batata tem uma variância muito superior às ### seguindo-se a variável milho, tendo as restantes variáveis ### variâncias quase nulas. ### a distância euclideana entre marcadores de indivíduos ### corresponde à distância de Mahalanobis entre os indivíduos. ### Como as distâncias de Mahalanobis são inferiores nas direções ### de maior variabilidade, as distâncias euclideanas entre os ### os marcadores de indivíduos na direção do PC1 ### são menores para refletir as distâncias de Mahalanobis ### ao longo do PC1, ocorrendo então uma compressão na direção eixo PC1 ### Finalmente a projeção ortogonal do marcador de cada individuo ### sobre a reta definida por uma marcador de variável é proporcional ### ao valor que essa variável toma nesse indivíduo ### (após divisão pelo desvio-padrão da variável). ### Conclui-se assim que os concelhos ### "Santarem", "Chamusca" , "Alpiarca", "V.N.Barquinha", ### "Coruche", "Entroncamento" e "Torres Novas" ### que se encontram do lado direito do biplot ### (identificando marcador da variável batata com o eixo PC1) ### têm uma produção por ha superior à média (12.28565), ### que é tanto maior quanto mais o concelho se encontrar à direita. ### Em sentido oposto, os concelhos do lado esquerdo do biplot têm uma ### produção por ha inferior à média, que é tanto menor quanto ### mais o concelho se encontrar para a esqueda... ### Uma análise semelhante pode ser feita em relação à produção ### de milho. Neste caso já não podemos identificar a reta definida ### pelo marcador da variável com o eixo PC2 e tem que se efetuar ### a projeção perpendicularmente à reta... head(X) colMeans(X) sort(X[,"batata"]) rownames(X)[order(X[,"batata"])] color_biplot(X,labels.ind=rownames(X),col.ind = as.numeric(X[,"batata"]>mean(X[,"batata"]))+3) color_biplot(X) sort(X[,"milho"]) rownames(X)[order(X[,"milho"])] color_biplot(X,labels.ind=rownames(X),col.ind = as.numeric(X[,"milho"]>mean(X[,"milho"]))+3) ### qualidade de representacao de cada concelho cos2=rowSums(scores[,1:2]^2)/rowSums(scores^2) plot(scores[,1:2],asp=TRUE,col="red",pch=16,cex=cos2) abline(h=0,lty=2,col="gray") abline(v=0,lty=2,col="gray") ####### ### b) View(X) View(Z) Z<-scale(X,scale=TRUE) Z_ACP.Z<-prcomp(X,scale=TRUE) ## PCA sobre matriz de correlaçoes ### PCA sobre dados normalizados = PCA sobre matriz de correlaçoes Z_ACP<-prcomp(Z,scale=FALSE) summary(Z_ACP) summary(X_ACP.Z) ### As primeiras explicam agora uma % muito inferior da variancia ### como é normal na ACP sobre uma matriz de correlação ### A núvem de pontos tende a ser mais arredondada para dados normalizados loadings.Z<-X_ACP.Z$rotation loadings.Z scores.Z<-X_ACP.Z$x scores.Z ## PC1 é essencialmente uma média ponderada de todas as produções ## explicando cerca de 1/3 da variabilidade total dos dados ## PC2 opõe os concelhos onde se produz mais fava, feijão e centeio, ## (Santarem, Rio Maior e Golegã) ## aos concelhos onde se produz mais trigo, aveia e cevada ## (Benavente, Coruche e Constancia) ## e portanto podemos pensar distinguindo cereais vs leguminosas plot(scores.Z[,1:2],asp=TRUE,col="red",pch=16) ## asp=TRUE -> keep the aspect ratio abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy text(scores.Z[,1],scores.Z[,2],labels=rownames(X),pos=1,cex=.5) ### nuvem de pontos é mais esférica que no caso ### da PCA sobre matriz de covariancias cos2.Z=rowSums(scores.Z[,1:2]^2)/rowSums(scores.Z^2) cos2.Z plot(scores.Z[,1:2],asp=TRUE,col="red",pch=16,cex=2*cos2.Z,main="qualidade de representação") abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy text(scores.Z[,1],scores.Z[,2],labels=rownames(X),pos=1,cex=.5) biplot(X_ACP.Z) abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy # Os marcadores de variaveis representados mais curtos estão mal # representados no PFP sort(X[,"fava"]) rownames(X)[order(X[,"fava"])] color_biplot(Z,labels.ind=rownames(X),col.ind = as.numeric((X[,"cevada"]>mean(X[,"cevada"]))&((X[,"grao"]>mean(X[,"grao"]))))+3) color_biplot(Z,labels.ind=rownames(X),col.ind = as.numeric(X[,"batata"]>mean(X[,"batata"]))+3) color_biplot(Z,labels.ind=rownames(X),col.ind = as.numeric(X[,"milho"]>mean(X[,"milho"]))+3) ## ii) cor(X,+scores.Z) loadings.Z cor(X,+scores.Z)/loadings.Z ### devolve para cada PC o sdev(PC)=sqrt(lambda) ### porquê? color_biplot(Z,labels.ind=rownames(X),col.ind = as.numeric(X[,"milho"]>mean(X[,"milho"]))+3,axes=c(3,6)) ### iii ) ######################################################## ### Exercício 10 ### a) summary(trigo) trigo.acpZ<-prcomp(trigo,scale=TRUE) summary(trigo.acpZ) ### alternativa equivalente trigo.Z<-scale(trigo,scale=TRUE) trigo.acpZ.bis<-prcomp(trigo.Z) loadings.trigo<-trigo.acpZ$rotation loadings.trigo scores.trigo<-trigo.acpZ$x scores.trigo ### b) cos2trigo=rowSums(scores.trigo[,1:2]^2)/rowSums(scores.trigo^2) plot(scores.trigo[,1:2],asp=TRUE,col="red",pch=16,cex=2*cos2trigo) ## asp=TRUE -> keep the aspect ratio abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy text(scores.trigo[,1],scores.trigo[,2],labels=rownames(trigo),pos=1) ## PC1 opõe os anos 1920-21 e 1927-28 ao ano 1929-30 ## vamos usar estas observações mais extremas ao longo da PC1 ## para caracterizar essa componente ## características dos anos x1, x2, x3,x4,x5 ## 1920-21 ## 1927-28 ## 1929-30 sort(trigo[,"x1"]) rownames(trigo)[order(trigo[,"x1"])] sort(trigo[,"x2"]) rownames(trigo)[order(trigo[,"x2"])] sort(trigo[,"x3"]) rownames(trigo)[order(trigo[,"x3"])] sort(trigo[,"x4"]) rownames(trigo)[order(trigo[,"x4"])] sort(trigo[,"x5"]) rownames(trigo)[order(trigo[,"x5"])] plot(scores.Z[,1:2],asp=TRUE,col="red",pch=16,cex=2*cos2.Z,main="qualidade de representação") abline(h=0,lty=2,col="gray") ## eixo dos xx a cinzento tracejado abline(v=0,lty=2,col="gray") ## idem para o eixo dos yy text(scores.Z[,1],scores.Z[,2],labels=rownames(X),pos=1,cex=.5) ### c) round(cor(trigo,scores.trigo),5) ### biplot(trigo.acpZ) color_biplot(trigo.Z,labels.ind=rownames(trigo),col.ind = as.numeric((trigo[,"x4"]>mean(trigo[,"x4"]))&(trigo[,"x5"]>mean(trigo[,"x5"]))&(trigo[,"x3"]mean(trigo[,"x3"])))+3)