# NOT RUN {
# Initial parameter values
beta <- c(-1,1.50); phi <- 5; rho <- 0.45; tau2 <- 0.80; sigma2 <- 1.5
# Simulating data
n1 <- 9 # Number of spatial locations
n2 <- 4 # Number of temporal index
set.seed(700)
x.coord <- round(runif(n1,0,10),9) # X coordinate
y.coord <- round(runif(n1,0,10),9) # Y coordinate
coordenadas <- cbind(x.coord,y.coord) # Cartesian coordinates without repetitions
coord2 <- cbind(rep(x.coord,each=n2),rep(y.coord,each=n2)) # Cartesian coordinates with repetitions
time <- as.matrix(seq(1,n2,1)) # Time index without repetitions
time2 <- as.matrix(rep(time,n1)) # Time index with repetitions
x1 <- rexp(n1*n2,2)
x2 <- rnorm(n1*n2,2,1)
x <- cbind(x1,x2)
media <- x%*%beta
# Covariance matrix
H <- as.matrix(dist(coordenadas)) # Spatial distances
Mt <- as.matrix(dist(time)) # Temporal distances
Cov <- CovarianceM(phi,rho,tau2,sigma2,distSpa=H,disTemp=Mt,kappa=0,type.S="exponential")
# Data
require(mvtnorm)
y <- as.vector(rmvnorm(1,mean=as.vector(media),sigma=Cov))
perc <- 0.2
aa=sort(y); bb=aa[1:(perc*n1*n2)]; cutof<-bb[perc*n1*n2]
cc=matrix(1,(n1*n2),1)*(y<=cutof)
y[cc==1] <- cutof
# Estimation
est_teste <- EstStempCens(y, x, cc, time2, coord2, inits.phi=3.5, inits.rho=0.5, inits.tau2=0.7,
type.Data="balanced", cens.type="left", method="nlminb", kappa=0,
type.S="exponential",
IMatrix=TRUE, lower.lim=c(0.01,-0.99,0.01), upper.lim=c(30,0.99,20), M=20,
perc=0.25, MaxIter=25, pc=0.2, error = 10^-6)
# }
# NOT RUN {
# New York data
library(spTimer)
library(sp)
library(rgdal)
# Transform coordinates
station <- as.data.frame(NYdata[,2:3])
names(station) <- c("lon","lat")
coordinates(station) <- ~lon+lat
proj4string(station) <- CRS("+init=epsg:32116 +proj=longlat +datum=NAD83 +no_defs
+ellps=GRS80 +towgs84=0,0,0")
station2 <- spTransform(station,CRS("+proj=utm +zone=18 +north +datum=NAD83"))
station <- as.data.frame(station2)
names(station) <- c("x.Coord","y.Coord")
NYdata <- cbind(NYdata,station)
coord <- unique(station)
EstEstimation <- c(1,2,3,5,6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,26,27,28)
EstEstimation <- coord[EstEstimation,]
dataEstimation <- NYdata[(NYdata[,11]%in%EstEstimation[,1])&(NYdata[,12]%in%EstEstimation[,2]),]
cc <- vector("numeric",length=nrow(dataEstimation))
cc[is.na(dataEstimation$o8hrmax)==1] <- 1
time <- rep(seq(1,62),nrow(EstEstimation))
dados <- cbind(dataEstimation,cc,time)
names(dados) <- c("s.index","Longitude","Latitude","Year","Month","Day","o8hrmax","cMAXTMP",
"WDSP","RH","x.Coord","y.Coord","censure","time")
# Data
y <- dados$o8hrmax
cc <- dados$censure
x <- as.matrix(cbind(dados$cMAXTMP,dados$WDSP,dados$RH))
tempo <- as.data.frame(dados$time)
coord <- as.matrix(cbind(dados$x.Coord/1000,dados$y.Coord/1000)) # Coordinates in Km
# Power exponential model for the spatial covariance structure
est <- EstStempCens(y, x, cc, tempo, coord, inits.phi = 50, inits.rho = 0.30, inits.tau2 = 25,
type.Data="balanced", cens.type="missing", method="nlminb", kappa=0.50,
type.S="pow.exp",
IMatrix=TRUE, lower.lim=c(0.01,-0.80,0.01), upper.lim=c(500,0.80,150), M=20,
perc=0.25, MaxIter=500, pc=0.2, error = 10^-5)
# }
Run the code above in your browser using DataCamp Workspace