### ### carregando os dados, disponivei em ### "http://www.leg.ufpr.br/doku.php/publications:papercompanions:pablfetal" dim(valencia <- read.table("valenciapab.txt", header=TRUE)) head(valencia) valencia$inplaaca <- ifelse(valencia$nacaplan>0,1,0) ### dimensão da lattice (ddd <- sapply(lapply(valencia[,1:3], unique), length)) ### organizando os dados selecionados me lattice (array) dim(irl <- tapply(valencia$inplsin, valencia[,1:3], as.integer)) ### definindo dados de novas incidencias (dim(dtmp <- data.frame(t(matrix(irl, ncol=ddd[3]))))) inew <- sapply(dtmp, function(x) { id <- which(x[1:(ddd[3]-1)]>0) x[id+1] <- NA x }) dim(inew) ### checando... dtmp[,844] inew[,844] table(as.integer(t(inew))[-(1:1160)], irl[,,2:23]) table(as.integer(t(inew))[-(1:1160)], irl[,,1:22]) ### carregando função de lag espaco temporal, disponivel em ### "http://www.leg.ufpr.br/doku.php/publications:papercompanions:pablfetal" source("../../R/slagfun.R") ### criando data.frame de dados para análise ### notar que C aqui é para vizinhas na ***mesma linha*** ### e R é para vizinhas na ***entre linha*** df <- data.frame(y=as.integer(t(inew)), snrl3(y~S+R+C,irl,-1), snrl3(y~S+R+C,irl,0), nacaplan=valencia$nacaplan, inplaaca=valencia$inplaaca) dim(df) names(df) ### checando table(df$y.S.0, df$inplaaca) table(valencia$inplsin, valencia$inplaaca) table(df[,1:2]) table(df[,c(1, 5)]) ### definindo observações a serem consideradas ### desconsiderando primeira avaliação como resposta e bordas id.df <- sort(intersect(which(!apply(is.na(df), 1, any)), array(1:prod(ddd), ddd)[2:57, 2:19, 2:23])) length(id.df) ### ajustando os modelos m <- list() m$m1 <- glm(y~y.C..1+y.R..1, binomial, df[id.df, ], model=T) m$m2 <- glm(y~y.C.0+y.R.0, binomial, df[id.df, ], model=T) m$m3 <- glm(y~inplaaca, binomial, df[id.df, ], model=T) m$m4 <- glm(y~inplaaca+y.C..1+y.R..1, binomial, df[id.df, ], model=T) m$m5 <- glm(y~inplaaca+y.C.0+y.R.0, binomial, df[id.df, ], model=T) m$m6 <- glm(y~nacaplan, binomial, df[id.df, ], model=T) m$m7 <- glm(y~nacaplan+y.C..1+y.R..1, binomial, df[id.df, ], model=T) m$m8 <- glm(y~nacaplan+y.C.0+y.R.0, binomial, df[id.df, ], model=T) sapply(m, function(x) length(resid(x))) ### checando novamente... names(m$m5$model) table(m$m5$model[,1:2]) table(valencia$inplsin[id.df], valencia$inplaaca[id.df]) ### resultados lapply(lapply(m, coef), round, 4) lapply(m, function(x) round(coef(summary(x))[,4],4)) round(sapply(m, AIC),2) sapply(m, logLik) cenarios <- expand.grid(incidac=0:1, mesmlinhaC=0:2, entrelinhaR=0:2) ### lembrando que R é **entre linha** e C é **mesma linha** coef(m$m5) eta.m5 <- apply(cenarios, 1, function(x) coef(m$m5)[1] + coef(m$m5)[2]*x[1] + coef(m$m5)[3]*x[2] + coef(m$m5)[4]*x[3]) coef(m$m1) eta.m1 <- apply(cenarios, 1, function(x) coef(m$m1)[1] + coef(m$m1)[2]*x[2] + coef(m$m1)[3]*x[3]) coef(m$m3) eta.m3 <- coef(m$m3)[1] + coef(m$m3)[2]*0:2 round(100/(1+exp(-tapply(eta.m5, cenarios[,3:1], mean))), 2) round(100/(1+exp(-tapply(eta.m1, cenarios[,-1], mean))), 2) round(100/(1+exp(-eta.m3)),2)