Hi,

I have a model, which has a logaritmic response, but there is an offset for this response.
I use optimize to find a value to minimise the residuals of the model:

f2 <- function(x,df) sum(summary(lm(log((y-x)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df))$residuals^2)
start <- optimize(f2,c(0,min(y)),df=df)$minimum
mod <- lm(log((y-start)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df)
As I have already some data, I have build my model, but I would like to simulate the effect of adding more data points on the confidence interval for the predicted values.

So my code was:

forcast <- function(mod2,start, df,stud.temp = 25,stud.end=52,sims=100,Temp=c(25,40), Mnd=c(0.5,1,1.5),sd=0.005){
    df <- as.data.frame(expand.grid(Mnd,Temp))
    names(df) <- c("Mnd","Temp")
    df$Pred <- exp(predict(mod2,list(Temp= df$Temp, Mnd=df$Mnd)))*df$Mnd+start
    res.int <- c(1:sims)
    f2 <- function(x,df) sum(summary(lm(log((y-x)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df))$residuals^2)
    for (i in 1:sims){
        y <- rnorm(nrow(df),df$Pred,sd)
        start <- optimize(f2,c(0,min(y)),df=df)$minimum
        mod <- lm(log((y-start)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df)
        pr <- (exp(predict(mod,list(Temp= stud.temp, Mnd=stud.end),interval="confidence"))*stud.end+start)
        res.int[i] <- (pr[,3]-pr[,2])/2
    }
    cat(Mnd,round(mean(res.int),3),"\n")
}

Where one can input the base model, which is used for the prediction of the values for the generated expand.grid dataframe.
For this predicted values, we use the rnorm function to simulate the error, and then the model is fitted for this values, and the width of the confidence interval at the given temperature and time point is calculated.

As I would like to investigate the effect of the number of months used on the width of the confidence interval, I would like to speed up the coding, so I can use more simulations to investigate the effect.


Thanks in advance



Bart

	[[alternative HTML version deleted]]

