#####################################
## Q11 — GLM, Seleção de Modelos ##
## Autora: Victória Souza ##
#####################################sp<-read.csv("dadosRiqAb.csv",header=T)Função: Ler o arquivo csv, “header=T” indica que a primeira linha contém os nomes das colunas, é o cabeçalho da base de dados. Criar o data frame sp.
str(sp)## 'data.frame': 21 obs. of 4 variables:
## $ coleta: chr "C1" "C1" "C1" "C2" ...
## $ lagoa : chr "Araruama" "Marica" "Saquarema" "Araruama" ...
## $ Riq : int 23 24 37 21 29 30 18 23 35 13 ...
## $ Ab : num 5042 38048 38739 9022 17509 ...Função: Exibir a estrutura de sp (número de linhas e colunas, fatores$, características dos dados (e.g. character (categórico), integer (número inteiro)))
names(sp)## [1] "coleta" "lagoa" "Riq" "Ab"Função: Listar os nomes das colunas do data frame sp.
attach(sp)Função: Tornar colunas de sp acessíveis pelo nome, sem a necessidade de especificar fatores (dados$).
amb<-read.csv("dadosAmb.csv",header=T)Função: Ler o arquivo csv, indicando que a primeira linha contém os nomes das colunas, é o cabeçalho da base de dados. Criar o data frame amb.
str(amb)## 'data.frame': 21 obs. of 9 variables:
## $ coleta: chr "C1" "C1" "C1" "C2" ...
## $ lagoa : chr "Araruama" "Marica" "Saquarema" "Araruama" ...
## $ subst : num 3.37 2.92 2.7 3.52 2.56 ...
## $ temp : num 25.9 25.1 24.2 27.7 29.1 ...
## $ ph : num 10.23 9.07 9.05 8.13 9.31 ...
## $ cond : num 790 355 562 784 410 ...
## $ sal : num 55 22.3 37.3 54.3 26.2 ...
## $ trans : num 100 50.3 70 53.7 53.8 ...
## $ prof : num 49.3 75 64.4 100 93.3 ...Função: Exibir a estrutura de amb (e.g. número de linhas e colunas, características dos dados (e.g. character, numeric))
names(amb)## [1] "coleta" "lagoa" "subst" "temp" "ph" "cond" "sal" "trans"
## [9] "prof"Função: Listar os nomes dos fatores (colunas) de amb.
amb2<-amb[-c(1,2,6)]Função: Criar o data frame amb2, removendo as colunas 1,2 e 6 de amb (fatores coleta, lagoa e cond)
names(amb2)## [1] "subst" "temp" "ph" "sal" "trans" "prof"Função: Verificar os nomes das variáveis restantes em amb2.
library(vegan)## Carregando pacotes exigidos: permuteFunção: Carregar o pacote vegan (padronização, distâncias, ordenação)
amb3<-decostand(amb2,method="standardize")Função: Estandardizar cada variável de amb2, gerando o data frame amb3. Evita que escalas diferentes dominem o ajuste.
library(car)## Carregando pacotes exigidos: carDatalibrary(lme4)## Carregando pacotes exigidos: MatrixFunção: Carregar os pacotes car (e.g. diagnóstico, vif, testes anova tipo II/III) e lme4 (e.g. modelos lineares, generalizados mistos).
m0<-glm(Riq~subst+temp+ph+sal+trans+prof, family=poisson, data=amb3) Função: Ajustar um modelo GLM (Modelo Linear Generalizado), utilizando família Poisson (ideal para dados de contagem), para riqueza (Riq) em função dos fatores fixos subst, temp, ph, sal, trans e prof, usando os preditores padronizadoos em amb3.
vif(m0)## subst temp ph sal trans prof
## 5.994489 1.241108 2.115595 5.964209 3.926874 1.712438Função: Calcular os fatores de inflação da variância (vif) para detectar colinearidade entre preditores em m0. (Regras de bolso: vif>5-10 indica colinearidade elevada, e o indicado é remover um dos preditores)
m1<-glm(Riq~temp+ph+sal+trans+prof, family=poisson, data=amb3) Função: ajustar um novo modelo GLM a partir do m0, removendo subst, para eliminar colinearidade e deixa-lo mais parcimonioso.
vif(m1)## temp ph sal trans prof
## 1.201130 1.222988 1.169516 1.642560 1.594999Função: reavaliar colinearidade no modelo m1.
m2<-glmer(Riq~temp+ph+sal+trans+prof+(1|lagoa), family=poisson, data=amb3)Função: ajustar um novo modelo, agora GLMM (Modelo Linear Generalizado Misto) a partir de m1. Agora com intercepto aleatório para lagoa, ou seja, considerando-o como um fator aleatório e controlando seu efeito.
anova(m1)## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Riq
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 20 39.337
## temp 1 0.0103 19 39.327 0.919160
## ph 1 1.6086 18 37.719 0.204694
## sal 1 5.1015 17 32.617 0.023906 *
## trans 1 7.4160 16 25.201 0.006465 **
## prof 1 7.6154 15 17.586 0.005787 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Função: Análise de Deviance Sequencial do GLM m1 para testar a contribuição de cada termo na ordem em que entra no modelo m1, por meio do teste qui-quadrado.
anova(m2)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## temp 1 1.3238 1.3238 1.3238
## ph 1 0.1668 0.1668 0.1668
## sal 1 0.2733 0.2733 0.2733
## trans 1 0.9913 0.9913 0.9913
## prof 1 8.8920 8.8920 8.8920Função: Exibir uma tabela de análise do GLMM ajustado, para os efeitos fixos.