Function to loop calculation of standard errors for model predictions-CodePudding

I am after some assistance to loop my code or create a function for my required calculations.

My data frame is as follows. All columns repeat the same value in every row except for newdat2$time, the value of which changes each row:

newdat2 <- data.frame(season = rep("Summer", 31), 
                      time = seq(0, 3, by = 0.1), 
                      temp = rep(21.79384, 31),
                      last.rain.bom = rep(4.232604, 31),
                      rain = rep(0.916501, 31),
                      wind = rep("nil", 31),
                      cloud = rep(40.20378, 31),
                      abundance = rep(117.6262, 31),
                      site = rep("Avalon", 31))

For each row of this data frame I would like to complete the below calculation. This calculation is calculating the standard errors of predictions from a fitted model, see here.

C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,time,21.8,4.23,0.917,0,0,0,40.2,4.78) # This represents covariate values of my fitted model. The value of time needs to change for each row of newdat2$time, all other values remain the same
s <- vcov(zib) # zib is my fitted model and this row of code is taking the variance covariance matrix of my fitted model. s is a matrix 27x27
newdat2$se <- sqrt(t(C) %*% s %*% C) # This then calculates the standard errors for my model predictions but C must change for each row of newdat2 to reflect the change in newdat2$time

For example, the first calculation completed by the loop/function would be

C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,0,21.8,4.23,0.917,0,0,0,40.2,4.78) # 0 is the first value of newdat2$time
s <- vcov(zib) 
newdat2$se <- sqrt(t(C) %*% s %*% C)

the second calculation completed by the loop/function would be

C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,0.1,21.8,4.23,0.917,0,0,0,40.2,4.78) # 0.1 is the second value of newdat2$time
s <- vcov(zib) 
newdat2$se <- sqrt(t(C) %*% s %*% C)

the third calculation completed by the loop/function would be

C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,0.2,21.8,4.23,0.917,0,0,0,40.2,4.78) # 0.2 is the third value of newdat2$time
s <- vcov(zib) 
newdat2$se <- sqrt(t(C) %*% s %*% C)

Any assistance to loop such a calculation or create a function that would enable this would be very much appreciated.

CodePudding user response：

I don't have the data or the expected result here, but this should work:
The idea is to make all the versions of the vector C into a matrix and then do the calculation with it. You would only need the diagonal elements of the resulting answer, so I think colSums(m * s %*% m) will give the same answer, but be faster.

C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,0,21.8,4.23,0.917,0,0,0,40.2,4.78)
m <- matrix(rep(C, length(newdat2$time)), ncol = length(newdat2$time))
m[19, ] <- newdat2$time
s <- vcov(zib)
newdat2$se <- sqrt(colSums(m * s %*% m))

This should be faster than looping.

CodePudding user response：

By looping, you can do it like the following:

newdat<-NULL
for(i in 1:length(newdat2$time))
{
    C = c(0,0,0,0,0,0,0.0,0,0,0,0,0,0, 0, 1,0,0,0,newdat2$time[i],21.8,4.23,0.917,0,0,0,40.2,4.78)
    s <- vcov(zib)
    newdat<-c(newdat,sqrt(t(C) %*% s %*% C))
}

Now you can just add the newdat vector to the dataframe. However, I agree with @Brian above that this one is slower as compared to the vectorised method he suggested.