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]]
