birth=matrix(data=NA,nrow=ngroup,ncol=nsample) for (g in 1:ngroup) ###compute the net births at each period {for (j in 1:nsample) { birth[g,j]=Nhat[g]*std_entry[g,j] } } B=vector(length=(ngroup*nsample)) for (g in 1:ngroup) {for (j in 1:nsample) {B[((g-1)*nsample+j)]=Nhat[g]*std_entry[g,j]} } ###B=vector(length=(ngroup*(nsample+1))) ###B[1:nsample]=birth[1,] ###B[(nsample+2):(ngroup*(nsample+1)-1)]=birth[2,] B b=vector(length=(ngroup*nsample)) for (g in 1:ngroup) {for (j in 1:nsample) {b[((g-1)*nsample+j)]=std_bstar[g,j]} } ###b[1:nsample]=std_bstar[1,] ###b[(nsample+2):(ngroup*(nsample+1)-1)]=std_bstar[2,] b partial=matrix(data=0,nrow=(ngroup*nsample),ncol=(nsample+1)) for (g in 1:ngroup) {for (i in ((g-1)*nsample+1):(g*nsample)) { partial[i,(nsample+1)]=B[i]/Nhat[g] } } for (i in 1:(nsample-1)) {partial[i,i]=B[i]/b[i]} if(GJSTL==1){ #Only for GJSTL model for (i in (nsample+1):(ngroup*nsample)) {for (j in 1:(nsample-1)) {if ((i-j)==nsample) {partial[i,j]=B[i]/b[i]} } } }#endif for (i in 1:nsample) {for (j in 1:nsample) {if (j>i) {partial[j,i]=-B[j]/(1-b[i])} } } if(GJSTL==1){ # Only for GJSTL model for (i in (nsample+1):(ngroup*nsample)) {for (j in (nsample+1):(ngroup*nsample)) {if (j>i) {partial[j,j-i]=-B[j]/(1-b[j-i])} } } }#endif H=matrix(data=0,nrow=(ngroup*nsample),ncol=nstd+1) H[,(nstd+1-(nsample+1)+1):(nstd+1)]=partial[,1:(nsample+1)] cov_net_birth=matrix(data=0,nrow=(ngroup*nsample),ncol=(ngroup*nsample)) for(g in 1:ngroup) { for (i in ((g-1)*nsample+1):(nsample*g)) {for (j in ((g-1)*nsample+1):(nsample*g)) {for (k in 1:(nstd+1)) {for (l in 1:(nstd+1)) {cov_net_birth[i,j]=cov_net_birth[i,j]+H[i,k]*H[j,l]*Cov_std1[k,l]} } } } } se_net_birth=vector(mode="numeric",length=(nsample*ngroup)) for(i in 1:(ngroup*nsample)) {se_net_birth[i]=sqrt(max(0,cov_net_birth[i,i]))} se_net_birth cov_pop_size=matrix(data=0,nrow=(ngroup*nsample),ncol=(ngroup*nsample)) se_pop_size=vector(mode="numeric",length=(nsample*ngroup)) H1=matrix(data=0,nrow=(ngroup*nsample),ncol=nstd+1) phi=vector(length=(ngroup*(nsample-1))) for (g in 1:ngroup) {for (j in 1:(nsample-1)) {phi[((g-1)*(nsample-1)+j)]=std_survival[g,j]} } phi pop_size=vector(length=(ngroup*nsample)) for (g in 1:ngroup) {for (j in 1:nsample) {pop_size[((g-1)*nsample+j)]=N_est[g,j]} } pop_size H1[1,]=H[1,] if(GJSTL==1){ H1[nsample+1,]=H[nsample+1,] }#endif for (i in 2:nsample) {H1[i,]=H[i,] for (j in 1:(nstd+1)) {H1[i,j]=H1[i,j]+H1[i-1,j]*phi[i-1] } H1[i,i-1]=H1[i,i-1]+pop_size[i-1] } if(GJSTL==1){ for (i in (nsample+2):(ngroup*nsample)) {H1[i,]=H[i,] for (j in 1:(nstd+1)) {H1[i,j]=H1[i,j]+H1[i-1,j]*phi[i-2] } H1[i,i-2]=H1[i,i-2]+pop_size[i-1] } }#endif for(g in 1:ngroup) { for (i in ((g-1)*nsample+1):(nsample*g)) {for (j in ((g-1)*nsample+1):(nsample*g)) {for (k in 1:(nstd+1)) {for (l in 1:(nstd+1)) {cov_pop_size[i,j]=cov_pop_size[i,j]+H1[i,k]*H1[j,l]*Cov_std1[k,l]} } } } } for(i in 1:(ngroup*nsample)) {se_pop_size[i]=sqrt(max(0,cov_pop_size[i,i]))} se_pop_size