In this exercise you will:
You need to have installed R, RStudio, and the necessary packages for
the course, including the ELMER
package. See how to
get set up for this course
'data.frame': 21 obs. of 2 variables:
$ Flow : int 24 36 24 25 72 37 46 24 37 46 ...
$ Problems: int 0 0 0 0 0 0 0 0 0 0 ...
A straight line did not model the mean of the river problems in the example found in ELMER at all well. Try fitting: \[\begin{aligned} \mu_{Y} & = & \beta_0 + \beta_1\ln{x} \\ var(Y) & = &k \mu_Y\end{aligned}\] Compare the fit of this model with the earlier one.
Step 1: Fit the unweighted model ModelU
, and make a copy
called ModelW
.
Step 2: Start to collect the values of the parameters you want to
track, and set up the Results
matrix to collate them.
Results<-NULL
Iteration=0
Beta0 <- coef(ModelW)[1]
Beta1 <- coef(ModelW)[2]
RMS<-sum(resid(ModelW)^2)/ModelW$df
Results<-rbind(Results,c(Iteration, Beta0, Beta1, RMS))
Step 3: Write out the working for the weighted model, and then wrap it in a loop to produce enough iterations. (10 should suffice).
You will need to:
update()
to incorporate the weights. for(i in 1:10){
Variance <-predict(ModelW)
Variance[Variance<=0]<-0.0001
ModelW<-update(ModelU, weight=1/Variance)
Iteration=i
Beta0 <- coef(ModelW)[1]
Beta1 <- coef(ModelW)[2]
RMS<-sum(1/Variance*resid(ModelW)^2)/ModelW$df
Results<-rbind(Results,c(Iteration, Beta0, Beta1, RMS))
}
Step 4: Print out the Results
to see that sufficient
iterations have ben completed to assure convergence.
(Intercept) log(Flow)
[1,] 0 5.9706779 -1.5518832 0.9502579
[2,] 1 0.1396734 -0.0335781 1.6046130
[3,] 2 3.2514315 -0.7631316 33.1513587
[4,] 3 3.3086476 -0.7737524 1.3057296
[5,] 4 3.3207104 -0.7765687 1.2674974
[6,] 5 3.3207242 -0.7765716 1.2628680
[7,] 6 3.3207257 -0.7765720 1.2628608
[8,] 7 3.3207257 -0.7765720 1.2628603
[9,] 8 3.3207257 -0.7765720 1.2628603
[10,] 9 3.3207257 -0.7765720 1.2628603
[11,] 10 3.3207257 -0.7765720 1.2628603
Step 5: Graph your new model.
RiverProb |> ggplot(aes(y=Problems,x=Flow)) + geom_point() +
xlab("Flow") + ylab("Number of problems") +
geom_function(fun=function(x) {Beta0 + Beta1*log(x)}, col="blue")