[R] AIC in MuMIn
Gavin Simpson
gavin.simpson at ucl.ac.uk
Tue Aug 17 15:27:34 CEST 2010
On Tue, 2010-08-17 at 16:05 +0800, elaine kuo wrote:
> Hello,
Why did you decide to post the exact same message from two different
email addresses??
> I am using package MuMIn to calculate AIC for a full model with 10
> explanatory variables.
> Thanks in advance in sharing your experience.
>
> Q1
> In the AIC list of all models, each model is differentiated by model number.
> Please kindly advise if it is possible to
> find the corresponding explanatory variable(s) for the model number.
Did you read ?dredge and look at the example to see how it worked? I've
never used this package and 30 seconds installing it, reading ?dredge
and running the very first example gave me the answer!
> # Example from Burnham and Anderson (2002), page 100:
> data(Cement)
> lm1 <- lm(y ~ ., data = Cement)
> dd <- dredge(lm1)
> subset(dd, delta < 4)
Global model: lm(formula = y ~ ., data = Cement)
---
Model selection table
(Int) X X1 X2 X3 X4 k R.sq Adj.R.sq RSS AIC AICc
11 52.58 1.468 0.6623 4 0.9787 0.9744 57.90 64.31 69.31
24 71.65 1.452 0.4161 -0.2365 5 0.9823 0.9764 47.97 63.87 72.44
23 48.19 1.696 0.6569 0.25 5 0.9823 0.9764 48.11 63.90 72.48
13 103.10 1.440 -0.6140 4 0.9725 0.9670 74.76 67.63 72.63
25 111.70 1.052 -0.41 -0.6428 5 0.9813 0.9750 50.84 64.62 73.19
delta weight
11 0.000 0.572
24 3.125 0.120
23 3.163 0.118
13 3.322 0.109
25 3.879 0.082
(or make the last line just read dd to print the full table)
So this clearly shows which variables (X, X1, etc) are in each model.
So, without anything further to go with regarding your problem, does
this answer this question?
> Q2 error message
> I tried to display sub-model with only temp_ran using the code below but
> failed.
> Please kindly suggest the potential failure cause.
>
> code
>
> rm(list=ls())
Please don't do the above; people might run your code and delete their
workspace if they aren't looking closely!
> library(MuMIn)
> datam <-read.csv("c:/migration/Mig_ratio_20100817.csv",header=T,
> row.names=1)
You do realise this isn't reproducible as I don't have your data file?
This is not what the posting guide asks you for. If you can't make the
problem appear with available data, then you might need to make *your*
data available for us to help.
> mig.stds <-glm(datam.std$SummerM_ratio~
> datam.std$temp_max++datam.std$evi_mean+datam.std$topo_var+datam.std$topo_mean+
> datam.std$coast+datam.std$Iso_index_0808,data=datam.std,family=gaussian)
Uggg! Ok, first lesson before we go any further; STOP using the data
frame name in your model formula. This is WRONG, WRONG, WRONG!!!! (Is
that clear enough ;-) I don't know where you learnt this, but unlearn
it, FAST! If someone taught you this, have a word with them.
Your model should be fitted like this:
mig.stds <-lm(SummerM_ratio ~ temp_max + evi_mean + topo_var +
topo_mean + coast + Iso_index_0808,
## now tell R were to find the variables in formula
data = datum.std)
## If you are fitting a Gaussian GLM it is better fitted with lm()
> target.model <- dredge(mig.stds, subset=datam.std$temp_ran)
But temp_ran is not in your model...
> error in eval(expr, envir, enclos), 'temp_ran' not found
Which would explain the error:
> dd2 <- dredge(lm1, subset = FOO)
Error in eval(expr, envir, enclos) : object 'FOO' not found
When used properly (none of this datam.std$ business), subset will do
what you want:
> dd2 <- dredge(lm1, subset = X1)
> dd2
Global model: lm(formula = y ~ ., data = Cement)
---
Model selection table
(Int) X X1 X2 X3 X4 k R.sq Adj.R.sq RSS
3 52.58 1.4680 0.6623 4 0.9787 0.9744 57.90
10 71.65 1.4520 0.4161 -0.2365 5 0.9823 0.9764 47.97
9 48.19 1.6960 0.6569 0.2500 5 0.9823 0.9764 48.11
5 103.10 1.4400 -0.6140 4 0.9725 0.9670 74.76
11 111.70 1.0520 -0.4100 -0.6428 5 0.9813 0.9750 50.84
6 52.53 0.15510 1.4670 0.6410 5 0.9798 0.9731 54.86
8 105.00 -0.16160 1.4380 -0.6379 5 0.9735 0.9647 71.91
15 62.41 1.5510 0.5102 0.1019 -0.1441 6 0.9824 0.9736 47.86
12 47.68 -0.05068 1.7240 0.6632 0.2801 6 0.9824 0.9735 47.93
13 70.97 0.01981 1.4520 0.4221 -0.2282 6 0.9823 0.9735 47.94
14 111.50 0.13500 0.9923 -0.4747 -0.6273 6 0.9818 0.9727 49.44
16 59.10 -0.01831 1.5930 0.5447 0.1453 -0.1124 7 0.9824 0.9698 47.85
2 71.93 1.51200 1.7300 4 0.6843 0.6212 857.40
1 81.48 1.8690 3 0.5339 0.4916 1266.00
4 72.35 2.3120 0.4945 4 0.5482 0.4578 1227.00
7 82.99 1.93300 1.0250 -0.7430 5 0.7048 0.6064 801.80
> Q3 showing sub-models for two assigned explanatory variables
> The manual only explains how to exclude two variables.
> Please advise how to contain submodels regarding certain two or more
> variables.
If I understand you (you want to select models that contain X1 and X2,
say), then you need to notice that the exclude example is just logical
indexing.
So
# exclude models with with both X1 and X2
dredge(lm1, subset = !X1 | !X2)
becomes
# include models with with both X1 and X2
dredge(lm1, subset = X1 & X2)
> dredge(lm1, subset = X1 & X2)
Global model: lm(formula = y ~ ., data = Cement)
---
Model selection table
(Int) X X1 X2 X3 X4 k R.sq Adj.R.sq RSS AIC
1 52.58 1.468 0.6623 4 0.9787 0.9744 57.90 64.31
4 71.65 1.452 0.4161 -0.2365 5 0.9823 0.9764 47.97 63.87
3 48.19 1.696 0.6569 0.2500 5 0.9823 0.9764 48.11 63.90
2 52.53 0.15510 1.467 0.6410 5 0.9798 0.9731 54.86 65.61
7 62.41 1.551 0.5102 0.1019 -0.1441 6 0.9824 0.9736 47.86 65.84
5 47.68 -0.05068 1.724 0.6632 0.2801 6 0.9824 0.9735 47.93 65.85
6 70.97 0.01981 1.452 0.4221 -0.2282 6 0.9823 0.9735 47.94 65.86
8 59.10 -0.01831 1.593 0.5447 0.1453 -0.1124 7 0.9824 0.9698 47.85 67.83
AICc delta weight
1 69.31 0.000 0.659
4 72.44 3.125 0.138
3 72.48 3.163 0.135
2 74.18 4.869 0.058
7 79.84 10.520 0.003
5 79.85 10.540 0.003
6 79.86 10.540 0.003
8 90.23 20.920 0.000
Notice how both the X1 and X2 columns always contain a coefficient. Or
for three variables:
> dredge(lm1, subset = X1 & X2 & X4)
Global model: lm(formula = y ~ ., data = Cement)
---
Model selection table
(Int) X X1 X2 X3 X4 k R.sq Adj.R.sq RSS AIC
1 71.65 1.452 0.4161 -0.2365 5 0.9823 0.9764 47.97 63.87
3 62.41 1.551 0.5102 0.1019 -0.1441 6 0.9824 0.9736 47.86 65.84
2 70.97 0.01981 1.452 0.4221 -0.2282 6 0.9823 0.9735 47.94 65.86
4 59.10 -0.01831 1.593 0.5447 0.1453 -0.1124 7 0.9824 0.9698 47.85 67.83
AICc delta weight
1 72.44 0.000 0.953
3 79.84 7.399 0.024
2 79.86 7.418 0.023
4 90.23 17.800 0.000
HTH
G
>
>
> Elaine
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Dr. Gavin Simpson [t] +44 (0)20 7679 0522
ECRC, UCL Geography, [f] +44 (0)20 7679 0565
Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/
UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
More information about the R-help
mailing list