set.seed(seed)
y0=rpois(Nt, lambda=N1*exp(-2.75+1.75*x+1*u))
# Negative control outcome in control region
zi=rpois(Nt, lambda=Ni*exp(-2.5*x-4+0.5*u ))
# Outcome of interest in control region
yi=rpois(Nt, lambda=Ni*exp(-3+0.75*x+0.75*u))
dat=data.frame(t=rep(1:Nt, 2), trt=rep(trt, 2), c=rep(c, 2), z=c(z1, zi), y=c(y1, yi), N=c(rep(500, Nt), rep(750, Nt)), reg=rep(1:2, each=Nt))
p1=with(dat[dat$reg==1 & dat$trt==1, ], sum(y)/sum(N))
p0=sum(y0[which(trt==1)])/(500*sum(trt))
return(list(rawdat=dat,
alpha=log(p1/p0),
ymod=summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef,
zmod=summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef ))
}
grbg=gen.dat(u=u, x=x)
names(grbg)
grbg$alpha
grbg$ymod
gen.dat=function(u, x, seed=8675309, N1=500, Ni=750){
trt=ifelse(u>2.1, 1, 0) # rbinom(Nt, 1, expit(-0.5 + c))
z1=rpois(Nt, lambda=N1*exp(2*x-4.75+0.95*u))
# Outcome of interest in intervention region
set.seed(seed)
y1=rpois(Nt, lambda=N1*exp(-2.75+1.75*x-.5*trt+1*u))
# make counterfactual
set.seed(seed)
y0=rpois(Nt, lambda=N1*exp(-2.75+1.75*x+1*u))
# Negative control outcome in control region
zi=rpois(Nt, lambda=Ni*exp(-2.5*x-4+0.5*u ))
# Outcome of interest in control region
yi=rpois(Nt, lambda=Ni*exp(-3+0.75*x+0.75*u))
dat=data.frame(t=rep(1:Nt, 2), trt=rep(trt, 2), c=rep(c, 2), z=c(z1, zi), y=c(y1, yi), N=c(rep(500, Nt), rep(750, Nt)), reg=rep(1:2, each=Nt))
p1=with(dat[dat$reg==1 & dat$trt==1, ], sum(y)/sum(N))
p0=sum(y0[which(trt==1)])/(500*sum(trt))
ymod=summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
zmod=summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
est.lfish=ymod["trt", "Estimate"]-zmod["trt", "Estimate"]
return(list(rawdat=dat,
alpha=log(p1/p0),
ymod=ymod,
zmod=zmod,
est.lfish=est.lfish))
}
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg$ymod
grbg$zmod
### Make a little function to do this more than once:
gen.dat=function(u, x, seed=8675309, N1=500, Ni=750){
trt=ifelse(u>2.1, 1, 0) # rbinom(Nt, 1, expit(-0.5 + c))
z1=rpois(Nt, lambda=N1*exp(2*x-4.75+0.95*u))
# Outcome of interest in intervention region
set.seed(seed)
y1=rpois(Nt, lambda=N1*exp(-2.75+1.75*x-.5*trt+1*u))
# make counterfactual
set.seed(seed)
y0=rpois(Nt, lambda=N1*exp(-2.75+1.75*x+1*u))
# Negative control outcome in control region
zi=rpois(Nt, lambda=Ni*exp(-2.5*x-4+0.5*u ))
# Outcome of interest in control region
yi=rpois(Nt, lambda=Ni*exp(-3+0.75*x+0.75*u))
dat=data.frame(t=rep(1:Nt, 2), trt=rep(trt, 2), c=rep(c, 2), z=c(z1, zi), y=c(y1, yi), N=c(rep(500, Nt), rep(750, Nt)), reg=rep(1:2, each=Nt))
p1=with(dat[dat$reg==1 & dat$trt==1, ], sum(y)/sum(N))
p0=sum(y0[which(trt==1)])/(500*sum(trt))
ymod=summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
zmod=summary(glm(z~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
est.lfish=ymod["trt", "Estimate"]-zmod["trt", "Estimate"]
return(list(rawdat=dat,
alpha=log(p1/p0),
ymod=ymod,
zmod=zmod,
est.lfish=est.lfish))
}
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
grbg=gen.dat(u=u, x=x)
grbg$alpha
grbg$est.lfish
mp=dat[dat$reg==1, ]
tmp$IR=with(dat[dat$reg==2, ], sum(y)/sum(N))
tmp$NIR=tmp$N*tmp$IR
summary(glm(y~offset(log(NIR)), data=tmp, family=poisson))$coef
tmp$JR=with(dat[dat$reg==2, ], sum(z)/sum(N))
summary(glm(y~offset(log(z*IR/JR)), data=tmp, family=poisson))$coef
log(0.5)
exp(-1.760516)
tmp=dat[dat$reg==1, ]
tmp$IR=with(dat[dat$reg==2, ], sum(y)/sum(N))
tmp$NIR=tmp$N*tmp$IR
summary(glm(y~offset(log(NIR)), data=tmp, family=poisson))$coef
tmp$JR=with(dat[dat$reg==2, ], sum(z)/sum(N))
summary(glm(y~offset(log(z*IR/JR)), data=tmp, family=poisson))$coef
exp(0.5)
plot(y1, ylim=range(c(z1, y1, yi, zi, y0)), pch=16, col="#00000075")
points(y0, col=ifelse(trt==1, "#ff0000", "#00000000"))
points(z1, col="#0000ff75", pch=16)
points(yi, ylim=range(c(zi, yi)), pch=1)
plot(y1, ylim=range(c(z1, y1, yi, zi, y0)), pch=16, col="#00000075", type="n")
points(yi, ylim=range(c(zi, yi)), pch=1)
points(zi, col="blue")
p1=with(dat[dat$reg==1 & dat$trt==1, ], sum(y)/sum(N))
p0=sum(y0[which(trt==1)])/(500*sum(trt))
p1/p0; log(p1/p0)
summary(glm(y~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
summary(glm(z~trt*as.factor(reg)+offset(log(N))-1, family='poisson', data=dat))$coef
library(sp)
?spplot
ni1=rbinom(1e3, 1, 0.2)
sum(ni1)
2/1e3
?rbinom
head(ni1)
table(ni1)
200/1e3
ni0=rbinom(1e3, 1, 0.1)
sum(ni1)
sum(ni0)
mean(rbinom(1e3, size=2e3, prob=0.3))
mean(rbinom(1e3, size=2e3, prob=0.3))
sum(c(rbinom(1e3, size=1, prob=0.2), rbinom(1e3, size=1, prob=0.1)))
sum(c(rbinom(1e3, size=1, prob=0.2), rbinom(1e3, size=1, prob=0.1)))
sum(c(rbinom(1e3, size=1, prob=0.2), rbinom(1e3, size=1, prob=0.1)))
sum(c(rbinom(1e3, size=1, prob=0.2), rbinom(1e3, size=1, prob=0.1)))
mean(rbinom(1e3, size=2e3, prob=0.15))
sum(c(rbinom(1.5e3, size=1, prob=0.2), rbinom(0.5e3, size=1, prob=0.1)))
(1.5*0.2+0.5*0.1)
(1.5*0.2+0.5*0.1)/2
mean(rbinom(1e3, size=2e3, prob=0.175))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
mean(rbinom(1e3, size=2e3, prob=(1.75*0.2+0.25*0.1)/2 ))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
sum(c(rbinom(1.75e3, size=1, prob=0.2), rbinom(0.25e3, size=1, prob=0.1)))
mean(rbinom(1e3, size=2e3, prob=(1.75*0.2+0.25*0.1)/2 ))
mean(rbinom(1e3, size=2e3, prob=(1.75*0.2+0.25*0.1)/2 ))
mean(rbinom(1e3, size=2e3, prob=(1.75*0.2+0.25*0.1)/2 ))
mean(rbinom(1e3, size=2e3, prob=(1.75*0.2+0.25*0.1)/2 ))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
sum(c(rbinom(1.75e4, size=1, prob=0.2), rbinom(0.25e4, size=1, prob=0.1)))
library(survey)
data(api)
dclus1<-svydesign(id=~dnum, data=apiclus1, fpc=~fpc)
dclus1
dim(dclus1)
table(apiclus1$fpc)
table(apiclus1$dnum)
length(apiclus1$dnum)
summary(dclus1)
dclus
dclus1
?svydesign
?api
names(apiclus2)
dclus2
dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
dclus2
table(apiclus2$dnum)
summary(dclus2)
svymean(~api00, dclus2)
svytotal(~enroll, dclus2, na.rm=TRUE)
table(apiclus2$fpc1)
table(apiclus2$fpc2)
svymean(~api00+meals+ell+enroll,dclus1,deff=TRUE)
svymean(~api00+meals+ell+enroll,dclus2,deff=TRUE,na.rm=TRUE)
dim(apiclus1)
dim(apiclus2)
?svymean
names(apiclus2)
dstrat<-svydesign(id=~1,strata=~stype,
#weights=~pw,
data=apistrat, fpc=~fpc)
svyratio(~api.stu,~enroll,dstrat)
sum(apipop$api.stu,na.rm=T)/sum(apipop$enroll,na.rm=T)
dstrat<-svydesign(id=~1,strata=~stype,
weights=~pw,
data=apistrat)#, fpc=~fpc)
svyratio(~api.stu,~enroll,dstrat)
dstrat<-svydesign(id=~1,strata=~stype,
weights=~pw,
data=apistrat, fpc=~fpc)
svyratio(~api.stu,~enroll,dstrat)
?srs
?sample
library(hhh)
library(surveillance)
?hhh4
data("influMen")
names(influMen)
head(influMen)
names(influMen)
names(influMen$observed)
names(influMen$head)
head(influMen$observed)
head(influMen$week)
head(influMen$state)
class(influMen)
dat=disProg2sts(influMen)
class(dat)
head(dat)
dat
names(dat)
names(dat@map)
dim(dat@map)
attributes(dat)
names(attributes(dat))
head(observed(dat))
head(state(dat))
head(dat@state)
install.packages(seegSDM)
install.packages("seegSDM")
install.packages("diseasemapping")
library(diseasemapping)
?kentucky
92/12
83636/2
83636/24
mu=2
alpha=-1
N=50
y1=rpois(1, N*exp(mu))
y2=rpois(1, N*exp(mu))
y3=rpois(1, N*exp(mu+alpha))
dat=data.frame(N=rep(N,3), y=c(y1, y2, y3), trt=c(0, 0, 1))
glm(y~trt+offset(log(N)), data=dat, family="poisson")
mod=glm(y~trt+offset(log(N)), data=dat, family="poisson")
summary(mod)
dat
(122/50)/((375+329)/(50+50))
log((122/50)/((375+329)/(50+50)))
log(0.15)
log(0.2)
0.15/0.20
Nis = 2 # number of areas
Ni = c(5000, 6000, 4000, 5500, 4500)[1:Nis]
runif(3, 4, 4)
runif(Nis, log(0.19), log(0.20))    # log Pr(ED | ILI, i)
epsi=runif(Nis, log(0.19), log(0.20))    # log Pr(ED | ILI, i)
deltai=runif(Nis, log(0.15), log(0.20))  # log Pr(ED | GI, i)
mu = -0.75 # log Pr(ILI)
nu = -0.5 # log Pr(GI)
alpha= -0.25 # log RR for treatment
dat=data.frame(i=rep(1:Nis, 2), Ni=rep(Ni, 2), c=rep(1:0, each=Nis), trt=rep(c(1, rep(0, Nis-1)), 2), muc=rep(c(mu, nu), each=Nis), epsci=c(epsi, deltai))
dat$trti=dat$trt*dat$c
dat$yc=with(dat, rpois(2*Nis, lambda=Ni*exp(epsci+muc+alpha*trti) ))
edif=c(with(dat[dat$c==1, ], epsci-epsci[1]), with(dat[dat$c==0, ], epsci-epsci[1]))
edif
dat
(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
=exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
edif
edif[(Nis+2):(2*Nis)]
Ni[2:Nis]
exp(0.23)
deltai=runif(Nis, log(0.195), log(0.205))  # log Pr(ED | GI, i)
dat=data.frame(i=rep(1:Nis, 2), Ni=rep(Ni, 2), c=rep(1:0, each=Nis), trt=rep(c(1, rep(0, Nis-1)), 2), muc=rep(c(mu, nu), each=Nis), epsci=c(epsi, deltai))
dat$trti=dat$trt*dat$c
dat$yc=with(dat, rpois(2*Nis, lambda=Ni*exp(epsci+muc+alpha*trti) ))
edif=c(with(dat[dat$c==1, ], epsci-epsci[1]), with(dat[dat$c==0, ], epsci-epsci[1]))
edif
(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
OI.factor=(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
NC.factor=exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
NC.factor/OI.factor
log(exp(alpha)*NC.factor/OI.factor)
log(NC.factor/OI.factor)
cat("Naive model bias: ", log(OI.factor))
cat("Adjusted model bias: ", log(NC.factor/OI.factor))
uniOI=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==1, ] )
summary(uniOI)
jtmod=glm(yc~offset(log(Ni))+c*trt, family=poisson, data=dat  )
summary(jtmod)
deltai=runif(Nis, log(0.19), log(0.20))  # log Pr(ED | GI, i)
dat=data.frame(i=rep(1:Nis, 2), Ni=rep(Ni, 2), c=rep(1:0, each=Nis), trt=rep(c(1, rep(0, Nis-1)), 2), muc=rep(c(mu, nu), each=Nis), epsci=c(epsi, deltai))
dat$trti=dat$trt*dat$c
dat$yc=with(dat, rpois(2*Nis, lambda=Ni*exp(epsci+muc+alpha*trti) ))
edif=c(with(dat[dat$c==1, ], epsci-epsci[1]), with(dat[dat$c==0, ], epsci-epsci[1]))
OI.factor=(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
NC.factor=exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
cat("Naive model bias: ", log(OI.factor))
cat("Adjusted model bias: ", log(NC.factor/OI.factor))
log(NC.factor/OI.factor)<log(OI.factor)
uniOI=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==1, ] )
summary(uniOI)
jtmod=glm(yc~offset(log(Ni))+c*trt, family=poisson, data=dat  )
summary(jtmod)
OI.factor
NC.factor
NC.factor/OI.factor
-0.1-0.25
-0.01-0.25
-0.07-0.25
uniNC=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==0, ] )
summary(uniNC)
epsi
deltai
epsi=runif(Nis, log(0.185), log(0.205))    # log Pr(ED | ILI, i)
dat=data.frame(i=rep(1:Nis, 2), Ni=rep(Ni, 2), c=rep(1:0, each=Nis), trt=rep(c(1, rep(0, Nis-1)), 2), muc=rep(c(mu, nu), each=Nis), epsci=c(epsi, deltai))
dat$trti=dat$trt*dat$c
dat$yc=with(dat, rpois(2*Nis, lambda=Ni*exp(epsci+muc+alpha*trti) ))
edif=c(with(dat[dat$c==1, ], epsci-epsci[1]), with(dat[dat$c==0, ], epsci-epsci[1]))
OI.factor=(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
NC.factor=exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
cat("Naive model bias: ", log(OI.factor))
cat("Adjusted model bias: ", log(NC.factor/OI.factor))
log(NC.factor/OI.factor)
uniOI=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==1, ] )
summary(uniOI)
uniNC=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==0, ] )
summary(uniNC)
jtmod=glm(yc~offset(log(Ni))+c*trt, family=poisson, data=dat  )
summary(jtmod)
OI.factor
NC.factor/OI.factor
1-NC.factor/OI.factor
1-OI.factor
epsi
exp(epsi)
exp(epsi)/exp(epsi[1])
1-exp(epsi)/exp(epsi[1])
crossprod(1-exp(epsi)/exp(epsi[1]), Ni)
crossprod(1-exp(epsi)/exp(epsi[1]), Ni)/Ni[2]
crossprod(1-exp(deltai)/exp(deltai[1]), Ni)/Ni[2]
dat=data.frame(i=rep(1:Nis, 2), Ni=rep(Ni, 2), c=rep(1:0, each=Nis), trt=rep(c(1, rep(0, Nis-1)), 2), muc=rep(c(mu, nu), each=Nis), epsci=c(epsi, epsi))
dat$trti=dat$trt*dat$c
dat$yc=with(dat, rpois(2*Nis, lambda=Ni*exp(epsci+muc+alpha*trti) ))
edif=c(with(dat[dat$c==1, ], epsci-epsci[1]), with(dat[dat$c==0, ], epsci-epsci[1]))
OI.factor=(exp(edif[2:Nis])%*%Ni[2:Nis])/sum(Ni[2:Nis])
NC.factor=exp(edif[(Nis+2):(2*Nis)])%*%Ni[2:Nis]/sum(Ni[2:Nis])
cat("Naive model bias: ", log(OI.factor))
cat("Adjusted model bias: ", log(NC.factor/OI.factor))
uniOI=glm(yc~offset(log(Ni))+trt, family=poisson, data=dat[dat$c==1, ] )
summary(uniOI)
jtmod=glm(yc~offset(log(Ni))+c*trt, family=poisson, data=dat  )
summary(jtmod)
dat
with(dat, yc/Ni)
dat$smr=with(dat, yc/Ni)
dat
with(dat, smr[1]/smr[2])
with(dat, log(smr[1]/smr[2]))
with(dat, log(smr[3]/smr[4]))
with(dat, log(smr[1]*smr[4]/(smr[2]*smr[3])))
with(dat, smr[3]/smr[4])
alpha
with(dat, yc[1]*yc[4]/(yc[2]*yc[3]))
with(dat, log(yc[1]*yc[4]/(yc[2]*yc[3])))
-0.2458124+0.08633648
setwd("~/Documents/Research/HFMD/BiometricsFinalPaper/WebSupplementary/SampleCode")
dat=read.table("SampleData.csv", header=T, sep=",")
alphaMat=do.call(rbind, by(dat, dat$strat, function(subdat){
tmp=data.frame(strat=subdat$strat[1])
tmp[, c("aES", "aCS", "aOS")]=colSums(subdat[, c("zES", "zCS", "zOS")])/colSums(dat[, c("zES", "zCS", "zOS")])
tmp[, c("aEM", "aCM", "aOM")]=colSums(subdat[, c("zEM", "zCM", "zOM")])/colSums(dat[, c("zEM", "zCM", "zOM")])
return(tmp)
}))
alphaMat$a.S=rowSums(alphaMat[, c("aES", "aCS", "aOS")])
alphaMat$a.M=rowSums(alphaMat[, c("aEM", "aCM", "aOM")])
dat$yES.hat=with(dat, zES+(yS-kS)*(alphaMat[dat$strat, "aES"]+zES)/(alphaMat[dat$strat, "a.S"]+kS) )
dat$yEM.hat=with(dat, zEM+(yM-kM)*(alphaMat[dat$strat, "aEM"]+zEM)/(alphaMat[dat$strat, "a.M"]+kM) )
dat$yCS.hat=with(dat, zCS+(yS-kS)*(alphaMat[dat$strat, "aCS"]+zCS)/(alphaMat[dat$strat, "a.S"]+kS) )
dat$yCM.hat=with(dat, zCM+(yM-kM)*(alphaMat[dat$strat, "aCM"]+zCM)/(alphaMat[dat$strat, "a.M"]+kM) )
dat$yOS.hat=with(dat, zOS+(yS-kS)*(alphaMat[dat$strat, "aOS"]+zOS)/(alphaMat[dat$strat, "a.S"]+kS) )
dat$yOM.hat=with(dat, zOM+(yM-kM)*(alphaMat[dat$strat, "aOM"]+zOM)/(alphaMat[dat$strat, "a.M"]+kM) )
dat$yE.hat=dat$yES.hat+dat$yEM.hat
dat$yC.hat=dat$yCS.hat+dat$yCM.hat
dat$yO.hat=dat$yOS.hat+dat$yOM.hat
dat$var.yES=with(dat, (yS-kS)*(alphaMat[dat$strat, "aES"]+zES)*(yS+alphaMat[dat$strat, "a.S"])*(alphaMat[dat$strat, "a.S"]+kS-alphaMat[dat$strat, "aES"]-zES)/((kS+alphaMat[dat$strat, "a.S"]+1)*(kS+alphaMat[dat$strat, "a.S"])^2) )
dat$var.yEM=with(dat, (yM-kM)*(alphaMat[dat$strat, "aEM"]+zEM)*(yM+alphaMat[dat$strat, "a.M"])*(alphaMat[dat$strat, "a.M"]+kM-alphaMat[dat$strat, "aEM"]-zEM)/((kM+alphaMat[dat$strat, "a.M"]+1)*(kM+alphaMat[dat$strat, "a.M"])^2) )
dat$var.yCS=with(dat, (yS-kS)*(alphaMat[dat$strat, "aCS"]+zCS)*(yS+alphaMat[dat$strat, "a.S"])*(alphaMat[dat$strat, "a.S"]+kS-alphaMat[dat$strat, "aCS"]-zCS)/((kS+alphaMat[dat$strat, "a.S"]+1)*(kS+alphaMat[dat$strat, "a.S"])^2) )
dat$var.yCM=with(dat, (yM-kM)*(alphaMat[dat$strat, "aCM"]+zCM)*(yM+alphaMat[dat$strat, "a.M"])*(alphaMat[dat$strat, "a.M"]+kM-alphaMat[dat$strat, "aCM"]-zCM)/((kM+alphaMat[dat$strat, "a.M"]+1)*(kM+alphaMat[dat$strat, "a.M"])^2) )
dat$var.yOS=with(dat, (yS-kS)*(alphaMat[dat$strat, "aOS"]+zOS)*(yS+alphaMat[dat$strat, "a.S"])*(alphaMat[dat$strat, "a.S"]+kS-alphaMat[dat$strat, "aOS"]-zOS)/((kS+alphaMat[dat$strat, "a.S"]+1)*(kS+alphaMat[dat$strat, "a.S"])^2) )
dat$var.yOM=with(dat, (yM-kM)*(alphaMat[dat$strat, "aOM"]+zOM)*(yM+alphaMat[dat$strat, "a.M"])*(alphaMat[dat$strat, "a.M"]+kM-alphaMat[dat$strat, "aOM"]-zOM)/((kM+alphaMat[dat$strat, "a.M"]+1)*(kM+alphaMat[dat$strat, "a.M"])^2) )
dat$var.yE=dat$var.yES+dat$var.yEM
dat$var.yC=dat$var.yCS+dat$var.yCM
dat$var.yO=dat$var.yOS+dat$var.yOM
dat$covECS.hat=with(dat, -(yS-kS)*(alphaMat[dat$strat, "aES"]+zES)*(alphaMat[dat$strat, "aCS"]+zCS)*(yS+alphaMat[dat$strat, "a.S"])/((kS+alphaMat[dat$strat, "a.S"]+1)*(kS+alphaMat[dat$strat, "a.S"])^2) )
dat$covECM.hat=with(dat, -(yM-kM)*(alphaMat[dat$strat, "aEM"]+zEM)*(alphaMat[dat$strat, "aCM"]+zCM)*(yM+alphaMat[dat$strat, "a.M"])/((kM+alphaMat[dat$strat, "a.M"]+1)*(kM+alphaMat[dat$strat, "a.M"])^2) )
pEj.hat=coef(glm(yE.hat~offset(log(Nj))-1+as.factor(strat), family=poisson, data=dat))
pCj.hat=coef(glm(yC.hat~offset(log(Nj))-1+as.factor(strat), family=poisson, data=dat))
pOj.hat=coef(glm(yO.hat~offset(log(Nj))-1+as.factor(strat), family=poisson, data=dat))
pEj.hat
with(dat[which(dat$strat==1), ], sum(yE.hat)/sum(Nj))
exp(pEj.hat)
yGt=0
yGt=data.frame(do.call(rbind, by(dat, dat$wk, function(wkdat){
c(wk=wkdat$wk[1], colSums(wkdat[, -c(1:2)]))
})))
yGt$EE.hat=rep(crossprod(exp(pEj.hat), Nj), Nt)
yGt$EC.hat=rep(crossprod(exp(pCj.hat), Nj), Nt)
yGt$EO.hat=rep(crossprod(exp(pOj.hat), Nj), Nt)
head(yGt)
dat[dat$wk==1, ]
Nj=dat[dat$wk==1, "Nj"]
yGt$EE.hat=rep(crossprod(exp(pEj.hat), Nj), Nt)
yGt$EC.hat=rep(crossprod(exp(pCj.hat), Nj), Nt)
yGt$EE.hat=rep(crossprod(exp(pEj.hat), Nj), dim(yGt)[1])
yGt$EC.hat=rep(crossprod(exp(pCj.hat), Nj), dim(yGt)[1])
yGt$EO.hat=rep(crossprod(exp(pOj.hat), Nj), dim(yGt)[1])
yGt$thetaEt.hat=with(yGt, yE.hat/EE.hat)
yGt$thetaCt.hat=with(yGt, yC.hat/EC.hat)
yGt$thetaOt.hat=with(yGt, yO.hat/EO.hat)
yGt$varEt=with(yGt, var.yE/(EE.hat^2))
yGt$varCt=with(yGt, var.yC/(EC.hat^2))
yGt$varOt=with(yGt, var.yO/(EO.hat^2))
yGt$covECt.hat=with(yGt, (covECS.hat+covECM.hat)/(EE.hat*EC.hat) )
yGt$lEt.hat=with(yGt, log(yE.hat/EE.hat))
yGt$lCt.hat=with(yGt, log(yC.hat/EC.hat))
yGt$lOt.hat=with(yGt, log(yO.hat/EO.hat))
yGt$var.lEt.hat=with(yGt, varEt/(thetaEt.hat^2))
yGt$var.lCt.hat=with(yGt, varCt/(thetaCt.hat^2))
yGt$var.lOt.hat=with(yGt, varOt/(thetaOt.hat^2))
yGt$covlECt.hat=with(yGt, covECt.hat/(thetaEt.hat*thetaCt.hat) )
library(INLA)
Nt=156 ## Number of times
temp=dat$temp[1:Nt] ## get temperature covariate
N <- 2*Nt
weights = rep(1, N)
data=list(lGt=c(yGt$lEt.hat, yGt$lCt.hat),
wkE=c(1:Nt, rep(NA, Nt)),
wkE2=c(1:Nt, rep(NA, Nt)),
wkC=c(rep(NA, Nt), 1:Nt),
wkC2=c(rep(NA, Nt), 1:Nt),
interE=c(rep(1,Nt), rep(NA,Nt)),
interC=c(rep(NA,Nt),rep(1,Nt)),
tempE=c(temp, rep(NA, Nt)),
tempC=c(rep(NA, Nt), temp)
)
for(j in 1:Nt) {
itmp = numeric(N)
itmp[] = NA
itmp[j] = 1
itmp[j+Nt] = 2
data = c(list(itmp), data)
names(data)[1] = paste("ii.", j, sep="")
}
add=""
for(j in 1:Nt) {
corr = yGt[j, "covlECt.hat"]/sqrt(prod(yGt[j, c("var.lEt.hat", "var.lCt.hat")]))
init.precE = log(1/yGt[j, "var.lEt.hat"])
init.precC = log(1/yGt[j, "var.lCt.hat"])
add = paste(add, paste(" +
f(", paste("ii.", j, sep=""), ", weights, model=\"iid2d\", n=2,
hyper = list(
prec1 = list(
initial =", init.precE,",
fixed = TRUE),
prec2 = list(
initial =", init.precC,",
fixed = TRUE),
cor = list(
initial = log((1+", corr, ")/(1-", corr, ")),
fixed = TRUE)))"))
}
S1=matrix(1:Nt, ncol=Nt, nrow=1)
formula.lGt = lGt ~ -1 + interE + interC + f(wkE, model="rw2", hyper=hypers[1], scale.model=F, constr=T, extraconstr=list(A=S1, e=0)) + f(wkC, model="rw2", hyper=hypers[2], scale.model=F, constr=T, extraconstr=list(A=S1, e=0)) + wkE2 + wkC2 + tempE + tempC
formula.lGt=update(formula.lGt, as.formula(paste(". ~ . ", add)))
mod= inla(formula.lGt, data=data,
family = "gaussian",
control.family = list(hyper = list(prec = list(initial = 10, fixed=T))),
control.predictor=list(compute=T),
control.compute=list(config=T),
control.inla=list(lincomb.derived.correlation.matrix=T),
control.fixed=list(prec=list(default=0.001), correlation.matrix=T) )
hypers.joint=c(list(prec = list(param = c(10, 1e-2) )),
list(prec = list(param = c(10, 1e-2) ))  )
mod= inla(formula.lGt, data=data,
family = "gaussian",
control.family = list(hyper = list(prec = list(initial = 10, fixed=T))),
control.predictor=list(compute=T),
control.compute=list(config=T),
control.inla=list(lincomb.derived.correlation.matrix=T),
control.fixed=list(prec=list(default=0.001), correlation.matrix=T) )
hypers=c(list(prec = list(param = c(10, 1e-2) )),
list(prec = list(param = c(10, 1e-2) ))  )
formula.lGt = lGt ~ -1 + interE + interC + f(wkE, model="rw2", hyper=hypers[1], scale.model=F, constr=T, extraconstr=list(A=S1, e=0)) + f(wkC, model="rw2", hyper=hypers[2], scale.model=F, constr=T, extraconstr=list(A=S1, e=0)) + wkE2 + wkC2 + tempE + tempC
mod= inla(formula.lGt, data=data,
family = "gaussian",
control.family = list(hyper = list(prec = list(initial = 10, fixed=T))),
control.predictor=list(compute=T),
control.compute=list(config=T),
control.inla=list(lincomb.derived.correlation.matrix=T),
control.fixed=list(prec=list(default=0.001), correlation.matrix=T) )
summary(mod)
