Computing std errors with R

[jlaake]

I posted a message in response to Carola's query and suggested an Excel sheet and showed the computation for RMark. You need not use RMark to take advantage of R to make these calculations. The following is R code for doing the calculation that is in the spreadsheet. Most of the code is assigning the values of the beta v-c matrix and beta values and outputting the results. The values can also input from a text file as shown in the second example that I constructed for Erin regarding her query a few months ago regarding average occupancy with a covariate. The lines starting with # are comments.

# R code for computing se(p*) for salamander example model
# assign values for beta v-c into a matrix
beta.vc=matrix(c(0.3369694, -0.305414194, -0.29150859, -0.29461909, -0.29971684,
-0.3054142, 0.688119206, 0.30309658, 0.30325514, 0.30351504,
-0.2915086, 0.303096577, 0.49276143, 0.30706109, 0.30518902,
-0.2946191, 0.303255138, 0.30706109, 0.49963140, 0.30481459,
-0.2997168, 0.303515038, 0.30518902, 0.30481459, 0.53105326),nrow=5,ncol=5,byrow=T)
# assign beta estimates for p to vector
beta=c(-1.5376878,-.3400082,1.1237206,0.9350627,0.5191252)
# create design matrix (standard treatment contrast) and compute p's using inverse
# logit link computed by function plogis
dm=matrix(c(rep(1,5), 0,1,0,0,0 ,0,0,1,0,0, 0,0,0,1,0, 0,0,0,0,1),nrow=5,ncol=5)
ps=plogis(dm%*%beta)
# compute derivative matrix of p's with respect to betas for logit link
D=dm*matrix(ps*(1-ps),nrow=length(ps),ncol=length(beta))
# compute v-c matrix of p's
V=D%*%beta.vc%*%t(D)
# compute p* = 1-(1-p1)*(1-p2)...
pstar=1-prod(1-ps)
# compute derivative vector of p* with respect to each p
dp=(1-pstar)/(1-ps)
# compute se of p*
sepstar=sqrt(t(dp)%*%V%*%dp)
# show values
pstar
sepstar


# Second example of computation of v-c matrix for Psi values where
# Psi is a function of site-specific covariates
# read in var-cov matrix for betas
beta_vc=as.matrix(read.table("beta_vc.txt",header=T))
row.names(beta_vc)=colnames(beta_vc)
# read in design matrix columns/rows for psi
dm=as.matrix(read.table("dm.txt",header=T))
# assign values for betas for psi
beta=c(3.117846,2.210670,-6.736828,11.743357,-2.520574)
# compute inverse logit link for psi estimates and show values
psi=plogis(dm%*%beta)
# compute derivative matrix of psi with respect to beta
D=dm*matrix(psi*(1-psi),nrow=length(psi),ncol=length(beta))
# compute v-c matrix for the psi parameters
vc.psi=D%*%beta_vc%*%t(D)
# compute std errors of psi which should match values in Presence output
se.psi=sqrt(diag(vc.psi))
# show table of psi estimates and std errors
print(data.frame(psi=psi,se=se.psi))
# show v-c matrix
print(vc.psi)
# compute var of mean prediction
cat("var of mean=",sum(vc.psi)/length(psi)^2,"\n")