|
##############Spatial Probit######################################################
rm(list=ls(all=TRUE))
campc<-read.table(“C:/Documents and
Settings//Desktop/ECVfamiglievaramb.txt”,col.names=c(“ID”,”xfam”,”yfam”,”Analfa”,”Istruzsup”,”Region”,
“Piso”,”acqua”,”servig”,”pareti”,”elettricit”,”basura”,”persc”,”TOTPER”,”TASAM”,”BEBES”,”Erosione”,
“Clima”,”Inondazioni”,”Vulcani”,”Areacod”,”REXPLRGE”,”Plinefamiglia”,”yprob”,”road1”,”road2”,”road3”,
“TOTAREA”,”Produttivit”,”PROT”,”IRR”,”Foreste”,”Arable”))
coord<-cbind(campc$xfam,campc$yfam)
coord<-coord/1000
cons<-log((campc$REXPLRGE))
n<-5630
media<-mean(cons)
dev<-var(cons)*n
diffm<-0
diffg<-0
k<-matrix(0,5630,1)
s<-matrix(0,5630,1)
s1<-0
ystar<-matrix(0,5630,1)
vic<-matrix(0,1,5630)
for (i in 1:5630) {for (j in 1:5630)
{dist<-sqrt((coord[i,1]-coord[j,1])^2+(coord[i,2]-coord[j,2])^2);
if (i!=j) {if (dist<=40) {vic[j]<-1/dist
diffm<-diffm+(((cons[i]-media)*(cons[j]-media))*vic[j]/dev)
diffg<-diffg+(((cons[i]-cons[j])^2)/dev)
k[j]<-k[j]+(1/dist)
s[i]<-s[i]+1/dist
s1<-s1+((1/dist)+(1/dist))^2}
else vic[j]<-0}
else vic[j]<-0}
ystar[i,1]<-vic%*%campc$yprob;
vic<-0}
a<-sum(k)
Moran<-(n/(a))*diffm
Geary<-(n-1)/(2*a)*diffg
s0<-sum(k)
s1<-(1/2)*s1
s2<-0
m4<-sum((cons-media)^4)/n
k1<-m4/((dev/n)^2)
for (i in 1:5630) {s2<-s2+(s[i]+k[i])^2}
varM<-((n*(((n)^2-3*n+3)*s1-n*s2+3*(s0)^2)-k1*(n*(n-1)*s1-2*n*s2+6*s0^2))/((n-1)*(n-2)*(n-3)*s0^2))-(1/(n-1)^2)
z<-(Moran-(-1/(n-1)))/sqrt(varM)
probabi<-dnorm(z)
ris<-c(“Moran=“,Moran)
if (probabi<0.05) ris1<-c(“Rifiuto H0: Not AUTOCORRELATION”)
else ris1<-c(“Accept H0: NON
AUTOCORRELAZIONE”)
ris2<-c(“Geary=“,Geary)
risultati<-cbind(ris,ris1,ris2)
write.table(risultati,file=“C:/Documents and Settings/Desktop/LMorant.txt”)
write.table(ystar,file=“C:/Documents and Settings/Desktop/tystar.txt”)
ystar<-matrix(scan(“c:/Documents and Settings/Desktop/tystar.txt”),ncol=1,byrow=T)
campc<-data.frame(campc,ystar)
campc$Erosione<-factor(campc$Erosione)
campc$Clima<-factor(campc$Clima)
campc$Arable<-factor(campc$Arable)
Spprobittot<-(glm(yprob~1+Analfa+Istruzsup+Piso+
acqua+servig+pareti+elettricit+basura+persc+TOTPER+TASAM+BEBES+Erosione+
Clima+Inondazioni+Vulcani+ystar+Areacod+road1+road2+road3+
TOTAREA+Produttivit+PROT+IRR+Foreste+Arable,data=campc,family=binomial(link=“probit”)))
deviance(Spprobittot)
summary(Spprobittot)
anova(Spprobittot,test=“Chisq”)
win.graph()
par(mfrow=c(2,2))
plot(resid(Spprobittot,type=“response”),xlab=“Y”,ylab=“ei”)
plot(resid(Spprobittot,type=“pearson”),xlab=“Y”,ylab=“epi”)
plot(resid(Spprobittot,type=“deviance”),xlab=“Y”,ylab=“edi”)
betapros<-c(Spprobittot$coef[1:21],Spprobittot$coef[23:33])
x<-matrix(scan(“c:/Documents and Settings/salvati/Desktop/xrural.txt”),ncol=32,byrow=T)
yprob<-x%*%betapros
povs<-pnorm(yprob)
plot(povs)
write.table(povs,file=“C:/Documents and Settings/salvati/Desktop/spatprobruramb.txt”)
########################Probit (Bigman)###########################################
rm(list=ls(all=TRUE))
campc<-read.table(“C:/Documents and
Settings/Desktop/ECVfamiglievaramb.txt”,col.names=c(“ID”,”xfam”,”yfam”,”Analfa”,”Istruzsup”,”Region”,
“Piso”,”acqua”,”servig”,”pareti”,”elettricit”,”basura”,”persc”,”TOTPER”,”TASAM”,”BEBES”,”Erosione”,
“Clima”,”Inondazioni”,”Vulcani”,”Areacod”,”REXPLRGE”,”Plinefamiglia”,”yprob”,”road1”,”road2”,”road3”,
“TOTAREA”,”Produttivit”,”PROT”,”IRR”,”Foreste”,”Arable”))
campc$Erosione<-factor(campc$Erosione)
campc$Clima<-factor(campc$Clima)
campc$Arable<-factor(campc$Arable)
probittot<-(glm(yprob~1+Analfa+Istruzsup+Piso+
acqua+servig+pareti+elettricit+basura+persc+TOTPER+TASAM+BEBES+Erosione+
Clima+Inondazioni+Vulcani+Areacod+road1+road2+road3+
TOTAREA+Produttivit+PROT+IRR+Foreste+Arable,data=campc,family=binomial(link=“probit”)))
deviance(probittot)
summary(probittot)
anova(probittot,test=“Chisq”)
win.graph()
par(mfrow=c(2,2))
plot(resid(probittot,type=“response”),xlab=“Y”,ylab=“ei”)
plot(resid(probittot,type=“pearson”),xlab=“Y”,ylab=“epi”)
plot(resid(probittot,type=“deviance”),xlab=“Y”,ylab=“edi”)
beta<-probittot$coef
x<-matrix(scan(“c:/Documents and Settings/Desktop/xurban.txt”),ncol=32,byrow=T)
yprob<-x%*%beta
povs<-pnorm(yprob)
plot(povs)
write.table(povs,file=“C:/Documents and Settings/Desktop/probruramb.txt”)
|
