######################################################################## ## Fri Sep 9 02:02:46 PDT 2011 ## This week we will: ## * gain greater facility with Emacs+R ## * gain experience with matrices and matrix multiplication ## * introduce ourselves to the "for() loop" ## * experience the wonders matplot() and barplot() ## * take a first look at apply() ## * discover a profound demographic truth ######################################################################## ## Comments on ## comments: ## While the main purpose of comments in your code is to help your ## biographer figure out how smart you were, they are also useful when ## you want to experiment by temporarily changing something -- in a ## way that you can change it back easily. For example: ## important.parameter<- 1 # old line of code important.parameter <- 2 # new line of code ## Only the second line will be executed by the R interpreter BUT if ## your program begins to emit smoke, you can easily restore it to it's ## previous state by 'uncommenting' the original line and deleting or ## commenting the new line ## For human readable comments, Recall that Emacs has ESC-Q -- which ## you can use to reformat or "fill" a bunch of raggedy comments. ## Emacs also has a "minor-mode" called "auto-fill-mode". You can ## toggle it on and off by typing C-x a. When auto-fill-mode is the ## word "Fill" will appear in the status bar-- and (more importantly) ## whenever you type past the end of a line, and Line-feed will be ## inserted along with the pound signs--IFF the current line begins ## with pound signs. It is *generally* useful, but can be annoying ## under some circumstances. Try toggling it on, and if it bugs you, ## toggle it off. ## One more nifty Emacs trick before we get to the R code: ## To highlight your copious comments with a whole line of ##, type ## M-72 # (that is the Meta key -- the the number 72 if you want 72 # ## signs, and then the characther that you want repeated. ######################################################################## ##### ## Now some R code ... ##### ######################################## ## the R commands below create a couple of vectors with very ## demographic sounding names ######################################## death.rate<-c(.0001,.01, .1, 1) birth.rate<-c(0, .75, .26, 0) init.pop<- c(1000,2000,1000,500) population<-array(NA,c(50,4)) # create an array called population, # fill it with NAs ## check the dimensionality of the population array dim(population) ## not too surprising ### how would you print population to the screen if you wanted to make ### sure that you really have created a 50X4 array of NAs ### Now let's do something moderately interesting with these objects: year<-1 ############################################# ##Question: What does this line do. If your not sure, how might you ##figure it out? ############################################# population[year,]<- init.pop ############################################# ## Question ## execute the next four lines and explain what you have accomplished ## by doing so ############################################# year<-year+1 population[year,1]<-sum(population[year-1,] * birth.rate) population[year,2:4]<- population[(year-1),1:3] * (1 - death.rate[1:3]) head(population) ## print the first 6 rows ############################################# ## Question ## compare the next four lines to the previous four lines ## what is the same about the code and what is different? ## ## What is the point of executing the next three lines? ## run the the next three lines anyway and examine the changed objects ## to see if you are right. ############################################# year<-year+1 population[year,1]<-sum(population[year-1,] * birth.rate) population[year,2:4]<- population[(year-1),1:3] * (1 - death.rate[1:3]) head(population) ## print the first 6 rows ############################################## ## Go ahead and run the previous 4 lines again ############################################# ############################################# ## Question ## What might the elements of the population array represent? ############################################# ######################################## ## There ought to be a way to do the sort of calcuclations ## we just did with out risking a repetative stress injury... ## and of course there is. We're going to cover loops and flow ## control in more detail later but the next few lines of code ## illustrate the "for loop" ## ######################################### ### the next 4 lines are some commands you have seen before. Here they ### are again to give you the bigger picture and make it easier ## to do the experiments called for later ## Repeated from above for completeness death.rate<-c(.0001,.01, .1, 1) # create death.rate vector birth.rate<-c(0, .75, .26, 0) init.pop<- c(1000,2000,1000,500) population<-array(NA,c(50,4)) # create an array called population, # fill it with NAs ### here is where the fun with loops begins population[1,]<- init.pop for (yindx in 2:50){ population[yindx,1]<-sum(population[yindx-1,] * birth.rate) population[yindx,2:4]<- population[(yindx-1),1:3] * (1 - death.rate[1:3]) } ######################################## ## Question ## Execute the commands above beginning with the one that asigns ## deat.rate all way to the } that closes the "for" loop. ## After you run those commands Print the population array ## What might the population array now contain? ##Specifically, what does the value in the 49'th row 3'rd column mean? ## HINT: Assume that the rate vectors hold age specific rates. ######################################## ############################################# ## Wouldn't it be nice to see these results graphically? ############################################# ## The plot() command is the simplest and most used plot in R, it ## draws scatter plots or line plots. plot(population[,2]) #draws a scatter plot of the second column (age #group ?) of the population matrix. By default, #the X values are just integers from #1:nrow(population). ## There are lot's of ways to improve the appearance of a plot() here ## are some: plot(population[,2],type='l', xlab="Period",ylab="Population in 2nd age group", main="Population forecast") footer<-paste("My Name",date()) # creates a text string which one # might tack on to the bottom of # graphs if one's name were "My # Name" ## mtext() aka "marginal text" puts text outside of the plot region. mtext(side=1,line=2,cex=.8,adj=1,text=footer) ## matplot() creates multiple scatter or line plots ######################################## ## Well, that was easy, but to understand what's going on, we'd like ## to see all of the age groups...wouldn't we? ######################################## ## Let's first choose some colors, one for each age group, most ## graphics functions will take a vector of colors pcolors<-c('slategray','darkorchid','chocolate3','#FF0000') #R #can undestand about 667 #colors by name, RGB #designations also work ## matplot is like executing plot() separately on each column of a ## matrix matplot(population,type='l',col=pcolors,xlab='Period',lwd=2,lty=1, ylab="Population by age group",main="Population forecast") ## create a legend legend(x='bottomright',fill=pcolors,legend=paste("age",1:4)) ############################################# ## Question ponder the graph you just drew, describe the population ## dynamics that it implies. ############################################# ## Let's look at a the result a few other ways... ## first, how about looking at total population. To get the total ## population at each period, we need to sum each row of the ## population matrix. apply() is the tool for this job total.pop<-apply(population,MARGIN=1,FUN=sum) length(total.pop) ## Perhaps the coolest feature of R graphics is that you can layer one ## thing on top of another. So assuming the graphic window still has ## the result of the most recent matplot(), we can use lines() to ## overlay the total population lines(total.pop,lty=2,lwd=3 ) ## R didn't complain... but it didn't seem to do anything either...? ## D'oh! min(total.pop) ## total.pop is off the chart. ## The easy way out is just to divide total.pop by 4 or so lines(total.pop/4,lty=2,lwd=2 ) ## But that would be the cowardly thing to do wouldn't it? ## A better alternative would be to cbind total.pop onto population ## and execute a new matplot() pcolors<-c(pcolors,'black') matplot(cbind(population,total.pop), type='l', col=pcolors,xlab='Period',lwd=2,lty=1, ylab="Population by age group",main="Population forecast") ## That's not all that satisfying either... How about a separate y ## axis on the right? ## the original matplot once again... pcolors<-c('slategray','darkorchid','chocolate3','#FF0000') matplot(population,type='l',col=pcolors,xlab='Period',lwd=2,lty=1, ylab="Population by age group",main="Population forecast") ##the par() function is used to adjust parameters that have to do with ## the behavior of the graphic window. par(new=TRUE) # counterintuitively tells the window NOT to clear # itself and redraw the axes when the next "high level" # plotting command is recieved. Without this command, # a plot() command would obliterate what was in the # graphics window and start again with new axes. plot(total.pop,type='l',lty=2,lwd=2,axes='n') #NOTE the axes='n' #argument ## Now to add a scale on the right side we use the axis() function axis(side=4) ## side 1=bottom, 2=left, 3=top, 4= right ## Hmmm this is misleading too isn't it. One might think that the ## rate of increase of total population exceeds the rate of increase ## of all of the age groups. Can that be? ############################################# ## QUESTION: create a plot somewhat like the one above using barplot() ## instead of matplot and showing the *proportion* of the total ## population that each age group comprises at each point in time ## rather than the levels. ## ## Note that barplot accepts a matrix just ## like matplot does but it creates a bar (or stacked bar) for each ## column. For this plot, you want 50 bars not 4, so you will need to ## use the t() function to transpose the population matrix. C-c C-v ## will "explain" barplot() ## If after a decent interval, you find this too frustrating to be ## useful, an answer is at the bottom of this file ## If on the other hand, you found the challenge exhilerating, then ## try superimposing on the barplot, a line showing the total ## population. par(mfrow=c(2,1)) # splits the graphics window horizontally (2 rowsX # 1 column) to allow for printing two graphs on one page ######################################## ## Question: You have no doubt noticed that when projecting the ## population in this way--with this particular initial population and ## these birth and death rates, in a few periods, the population ## settles into a constant rate of growth and a fixed age ## distribution. ## WITHOUT CHANGING EITHER THE BIRTH.RATES OR DEATH.RATES VECTORS, can ## you find a NONZERO starting population (init.pop) that causes the ## population to die out or explode in 50 periods? or to equilibrate ## with a different age break down. In other words, with fertility and ## mortality rates constant can you find a starting population (with ## NO zero elements) that will make the population increase at an ## increasing rate (explode) or die out? ## In order to experiment, you will want to create a new buffer/file ## that contains all the code necessary to run a trial with a set of ## init.pop that sets the initial population, and the birth and death ## rates and then runs the for() loop that populates the population ## matrix and then draws a graph. You can then execute all the lines ## in the buffer by typine C-c C-b; observe the graph that results; ## change init.pop (commenting out previous values..); rinse and ## repeat. ############################################# ## Another way to run this simulation is with matrix multiplication ## Mathematicians are fond of matrix algebra because one symbol ## replaces a whole mess of tedious calculations. If you don't ## remember your Linear Algebra, Wikipedia can refresh your memory on ## how it works. ############################################# Leslie<-rbind(birth.rate,diag(1-death.rate))[-5,] Leslie Leslie %*% init.pop head(population) ## Consider (Leslie%*% Leslie) %*% init.pop (Leslie%*% Leslie %*% Leslie) %*% init.pop ( (Leslie%*% Leslie %*% Leslie) %*% Leslie) %*% init.pop head(population) ## Question which row of population should equal this: (Leslie^10) %*% init.pop ## Question write a for loop that will produce Leslie raised to the ## 10th power -- in a matrix multiplication sense. Test your result ## agains the 11th row of population ## For the profoundly frustrated, an answer is at the bottom of this file ## That's it. Save your work -- but NOT your R workspace. ######################################################################## ## Frustration Management ######################################################################## ## This will create a barplot that shows the proportion breakdown of ## the simulated population across time. ## barplot(t(population/apply(population,1,sum)),col=pcolors,xlab='Period', ## ylab='proportion of population',ylim=c(0,1.2)) ## legend(x="top",legend=paste("age",1:4),ncol=4,bty='n', ## fill=pcolors) ## mtext(side=1,line=2,cex=.8,adj=1,text=footer) ## Raising square matrices to powers ## answer<-Leslie ## power<-10 ## for(i in 1:(power-1)){ ## answer<-answer%*%Leslie ## print(answer) ## } ## answer%*%init.pop