## Script Aula 13 - Aedes - buffers de 30 m #getwd() #setwd("C:/f_chiara_neto_ativos/cursos/aulas_FSP/a_pos/2019/PSP5124/aulas/aula11") load("aedes_acp_30m.Rdata") # INLA ## instalacao do INLA pela orientacao do livro Blangiardo e Cameletti #source("http://www.math.ntnu.no/inla/givemeINLA.R") ### instalando ultima versao mais estavel (a anterior a ultima) do pacote INLA #install.packages("INLA", repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE) ## instalando a versao ainda em teste (a ultima) #install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/testing"), dep=TRUE) #inla.upgrade(testing = F) # atualizar o pacote INLA #remove.packages("INLA") # remover o pacote INLA #carregando pacotes library(INLA) library(lattice) # para funcao 'dotplot' source("HighstatLibV10.R") # DATABASE banco <- read.csv("./bancos_aula13/banco_30m_acp.csv",header=T,sep=";") head(banco) str(banco) ############################################################################################# #### Resultados das Analise de Componentes Principais ### feita com o pacote 'psycho' # Principal Components Analysis # Call: principal(r = banco.pdr[, -c(1:11, 16, 17, 18, 19, 22)], nfactors = 3, # rotate = "varimax") # Standardized loadings (pattern matrix) based upon correlation matrix # RC1 RC2 RC3 h2 u2 com # AGUA 0.67 0.01 0.26 0.51 0.49 1.3 # ARVORE 0.80 0.25 0.03 0.71 0.29 1.2 # ASFALTO -0.82 0.22 0.10 0.73 0.27 1.2 # CIMENTO -0.35 0.64 0.46 0.75 0.25 2.4 # SOMBRA 0.46 0.33 0.60 0.68 0.32 2.5 # SOLO_EXPOSTO 0.01 -0.80 0.12 0.66 0.34 1.0 # TELHA_TELHA -0.22 -0.80 -0.12 0.70 0.30 1.2 # telha_am_log -0.41 0.59 -0.42 0.69 0.31 2.6 # grama_log 0.87 -0.09 0.02 0.76 0.24 1.0 # laje_log 0.02 0.11 -0.84 0.73 0.27 1.0 # # RC1 RC2 RC3 # SS loadings 3.06 2.30 1.56 # Proportion Var 0.31 0.23 0.16 # Cumulative Var 0.31 0.54 0.69 # Proportion Explained 0.44 0.33 0.23 # Cumulative Proportion 0.44 0.77 1.00 ############################################################## ### exploratory analysis ##### outliers MyVar <- c("RC1","RC2","RC3") Mydotplot(banco[,MyVar]) # usando a funcao 'HighstatLibV9.R' ### colinearidade - junto com os modelos com covar MyVar <- c("RC1","RC2","RC3") Mypairs(banco[,MyVar]) # sao componentes principais - e sao ortogonais ### relacao X x Y MyX <- c("RC1","RC2","RC3") Myxyplot(banco, MyX, "num_aeg_f3", MyYlab = "num aeg femeas") ################################################################### ###### standardizing covariates - resultados na escala de desvio padroes ### aqui nao eh necessario pq as componentes j? est?o padronizadas # str(banco) # head(banco) # banco.pdr <- banco # banco.pdr[,c(6:8)] <- scale(banco[,c(6:8)],scale =T) # str(banco.pdr) # head(banco.pdr) ######################### #### MODELOS str(banco) Loc <- cbind(banco$LONG_UTM, banco$LAT_UTM) Loc # dist entre os pontos D <- dist(Loc) par(mfrow = c(1,2), mar = c(5,5,2,2), cex.lab = 1.5) hist(D, freq = TRUE, main = "", xlab = "Distance between sites (km)", ylab = "Frequency") text(80, 470, "A", cex = 1.5) plot(x = sort(D), y = (1:length(D))/length(D), type = "l", xlab = "Distance between sites (km)", ylab = "Cumulative proportion") text(10, 1, "B", cex = 1.5) # mesh mesh <- inla.mesh.2d(Loc, max.edge=70, offset=100, cutoff=0) mesh$n # 984 par(mfrow = c(1,1)) plot(mesh) points(Loc, pch=8, cex=1.5, col=2) # Projector matrix A2 <- inla.spde.make.A(mesh, loc = Loc) dim(A2) ##### spde with PC priors spde.pc.priors <- inla.spde2.pcmatern(mesh, alpha = 2, prior.range=c(100, 0.5), prior.sigma=c(1, 0.5), constr=TRUE) w.index <- inla.spde.make.index( name = 'w', n.spde = spde.pc.priors$n.spde, n.group = 1, n.repl = 1) # Create a data frame with an intercept and covariates. N <- nrow(banco) ##################################################### ### 1.1 model only with the intercept and the spatial componente W str(banco) X <- data.frame(Intercept = rep(1, N)) #Stack Stk2 <- inla.stack( tag = "Fit", data = list(y = banco$num_aeg_f3), A = list(A2, 1), effects = list( w = w.index, #Spatial field X = as.data.frame(X))) #Covariates f.inter <- y ~ -1 + Intercept + f(w, model=spde.pc.priors) I.inter <- inla(f.inter, family = "poisson", data = inla.stack.data(Stk2), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stk2))) I.inter$dic$dic # 188.71 summary(I.inter) ####################### ################################################################## '' ### Modelos as componentes prinicipais - RC1, 2 and 3 X.RC1.2.3 <- data.frame(Intercept = rep(1, N), RC1 = banco$RC1, RC2 = banco$RC2, RC3 = banco$RC3) #Stack Stk.RC1.2.3 <- inla.stack( tag = "Fit", data = list(y = banco$num_aeg_f3), A = list(A2, 1), effects = list( w = w.index, #Spatial field X.RC1.2.3 = as.data.frame(X.RC1.2.3))) #Covariates f.RC1.2.3 <- y ~ -1 + Intercept + RC1 + RC2 + RC3 + f(w, model=spde.pc.priors) I.RC1.2.3 <- inla(f.RC1.2.3, family = "poisson", data = inla.stack.data(Stk.RC1.2.3), control.compute = list(dic = TRUE), control.predictor = list(A = inla.stack.A(Stk.RC1.2.3))) I.RC1.2.3$dic$dic # 187.37 summary(I.RC1.2.3) ### Betas # Fixed effects: # mean sd 0.025quant 0.5quant 0.975quant mode kld # Intercept 0.015 0.180 -0.361 0.023 0.347 0.039 0 # RC1 -0.068 0.156 -0.370 -0.069 0.245 -0.073 0 # RC2 -0.088 0.159 -0.403 -0.087 0.222 -0.086 0 # RC3 -0.549 0.138 -0.819 -0.550 -0.274 -0.552 0 ##################################################################### ### salvando a workspace save.image("aedes_acp_30m.Rdata") ################################################################