title: Which values to use for scaling out of sample PCA data

text: I have centered and scaled inputs via `prcomp ()`:               prOut<-prcomp(trainSet[,2:4],scale = TRUE,scores=TRUE)      I now want to use my completed model on new (future) data. I assume the correct approach is to use the training data `prOut$scale` and `prOut$center` values and apply them to prior to calculating the principal component scores for my new data? It doesn't seem right to scale the new data using the scale and center values for the new data set. thanks for any insight you can lend.


title: Linear Discrimination Scores

text: Suppose we have 4 groups of animals (pig, dog, rabbit, cat). There are 10 animals in each group. We measure the weight and pulse of each animal. ![enter image description here](http://i.stack.imgur.com/0ramF.jpg) Consider the following procedure in SAS:               proc discrim data=animals pool=test crossvalidate testdata=animals testout=a;     class group;     var weight pulse;     run;      Why wouldn't linear discriminant functions be returned in SAS?


title: Imprecise measurements and observation model

text: I am using sequential filtering framework to estimate the average speeds on different freeway segments from spot speeds obtained from multiple vehicles. Consider a mile-long stretch of a freeway, and a 2-minutes interval. I get speed measurements from $n$ vehicles (from an app like Waze) on this segment during this time interval. I estimate the average speed on the segment as the mean of these $n$ measurements. Let us call this estimate $\bar{v}$. I, somehow, in retrospect know the actual average speed $v$ on the segment. So, I defined an observation model as $$ \bar{v} = \alpha v + \beta + \mathcal{N}(0,\sigma^2)$$ Here $\alpha$ and $\beta$ are scalars to be determined, as is $\sigma$. This is just a linear regression model between $\bar{v}$ and $v$ and so I can estimate these parameters. No problem so far! BTW, the way my system is implemented, I **have** to find $\bar{v}$ and then use it as an observation to estimate $v$. I cannot directly use the $n$ speed measurements to estimate $v$. Now, I know that $\bar{v}$ was computed from $n$ measurements. So it matters whether $n = 5$ or $n = 30$. I would have more confidence in the value of $\bar{v}$ obtained from 30 measurements than just 5 measurements. I want to use this information in my observation model. In other words, I want to now treat $\bar{v}$ as an _imprecise_ measurement. So instead of a single value for $\bar{v}$, I get the value as well as a standard deviation $\delta$ (hence assuming $\bar{v}$ is normally distributed with variance $\delta^2$). $\delta$ value is known. How can I incorporate this information in the observation model and how do I now compute $\alpha$, $\beta$ and $\sigma$?


title: Variance of EM mean estimates in a simple mixture of two normals

text: Consider a mixture of two normal distributions: $ f(x) = p N(x|u_1, S_1) + (1-p) N(x|u_2, S_2) $ where N() is the normal pdf. $p$, $S_2$, and $S_2$ are known. The means are not. You can get the MLE of $u_1$ and $u_2$ from the EM algorithm. My question is, for a sample of size n, what are the variances of the EM estimates? My intuition tells me $S_1/(n \times p)$ and $S_2/(n \times (1-p))$ respectively, but I am not sure. Further, do the two estimates have a covariance? What would it be?


title: Dividing/Multiplying Gaussians

text: One of the steps in Expectation Propagation based inferencing is to be able to divide out a factor which has a normal form from the current posterior distribution (which is a high dimensional Gaussian). The idea is that the division removes the influence of that factor from the posterior distribution. So, assume that the posterior is multivariate normal with a n-dimensional mean and a $n$ $\times$ $n$ covariance matrix. Now, my likelihood depends on a subset of this n-dimensional parameters (actually 3. So, displacement along the 3 spatial dimensions and has a variance of $\sigma$). Just to make it clear, the likelihood is defined in terms of the residual as: $$ L(e) \propto \exp^{-0.5 e \sigma e} $$ where $e$ is the residual that depends non-linearly on sub-set of posterior parameters $w$ and has a precision given by $\sigma$ ($\sigma$ is a global noise precision). Just to explain, regardless of the fact that the likelihood is normally distributed wrt to the residuals it should still be seen as a function of $w$ for EP to work. Now, my question is how can I generate such a factor? So, the factor I need to be designed should be a normal distribution that is a function of all of $w$ but depends on only a subset of $w$ (those 3 elements). Also, I need to somehow generate a precision matrix that will be equivalent to the $\sigma$ and more importantly if I divide the posterior distribution by this Gaussian factor, it would be equivalent to removing the influence of this factor from the posterior distribution. I hope I have managed to explain the problem clearly. Please let me know if I can make it clearer.


title: Variance of arrival process with shifted exponential distribution

text: Here we have a arrival process. The inter-arrival time follows a shifted negative exponential distribution. The density function of the distribution is: $$f(t)=\lambda e^{-\lambda(t-\theta)},\quad\text{ where }t\ge\theta$$ How to derive the variance of the number of arrivals in time period of $T$?


title: What is the correct way to identify overdispersion?

text: I calculated the dispersion parameters for each of two blocks in a hierarchical logistic regression. I am not quite sure how to interpret the dispersion parameter. First I calculated the parameters by taking the chi-square goodness-of-fit statistic (hosmer-lemeshow-test) and dividing it by the degrees of freedom. For both blocks the values were less than 1. In the second step (to support the prior calculation) I took the 2-log- likelihood value and divided it by the degrees of freedom. In this case the parameters for both blocks were above 10. So is the first or the second calculation right? What can I do if there overdispersion?


title: Tobit or OLS question

text: My dependent variable is mostly continuous and positive, but has a modest number of zeros (10% of the sample). The results from tobit and ols are very similar. How can I formally compare the tobit and ols results? I think that I should be using ols, because tobit doesn't change anything.


title: Combine multiple boxplots in a lattice

text: I have a data-table that has about 26000 rows and about 35 columns. The columns are paired, so the values in columns 6 and 7 (for example) are related to each other, so are 8 and 9 and so on. There are 23 different types of annotations in the table, which I have read in as "factor". The ratio of these pairs of columns gives me a meaningful number, that I have to plot for each of the annotation. I was wondering if there is any way to have a lattice plot that will have say 15 boxplots in each panel, and 23 panels one for each annotation? UPDATE: Sample table.               structure(list(chromosome = structure(c(1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L), .Label = c("chr1", "chr2", "chr3"), class = "factor"),          start = c(1, 1, 1, 5663, 5726, 6360, 7548, 7619, 11027, 12158         ), end = c(5662, 7265, 5579133, 7265, 6331, 6755, 12710,          9274, 11556, 12994), strand = structure(c(1L, 1L, 3L, 1L,          1L, 1L, 3L, 3L, 1L, 3L), .Label = c("-", ".", "+"), class = "factor"),          annotation = structure(c(4L, 13L, 8L, 2L, 13L, 18L, 18L,          13L, 12L, 13L), .Label = c("3'-UTR", "5'-UTR", "BLASTN_HIT",          "CDS", "CDS_motif", "CDS_parts", "conflict", "Contig", "intron",          "LTR", "misc_feature", "misc_RNA", "mRNA", "polyA_site",          "promoter", "real_mRNA", "rep_origin", "repeat_region", "repeat_unit",          "rRNA", "snoRNA", "snRNA", "tRNA"), class = "factor"), Abp1D.sense = c(274.043090077,          222.027002967, 273.083037487, 38.3559401569, 80.7384755736,          15.9496926371, 54.9087080745, 127.744117176, 11.7165833969,          96.1925577965), Abp1D.antisense = c(125.681512904, 151.232091139,          254.813202986, 241.034453038, 84.3769908653, 199.467664241,          54.1912835565, 94.2017362521, 66.5142677515, 63.28607875),          Iki3D.sense = c(1214.1686727, 969.99693773, 261.416187303,          107.770848316, 151.518863438, 55.9449713698, 66.0800496533,          144.470307921, 21.9708783825, 52.6163190329), Iki3D.antisense = c(786.364743311,          728.647444388, 248.288893165, 523.636519401, 263.419180997,          351.558399018, 73.754086788, 130.973198864, 93.7873464478,          30.858803946), Iki3D.Rrp6D.sense = c(3068.90441567, 2486.4012139,          278.274812147, 428.928792511, 639.682546716, 134.968168726,          223.376134645, 491.4747595, 72.255001742, 201.429779476),          Iki3D.Rrp6D.antisense = c(1928.37423684, 1764.06364622, 271.050084744,          1181.76403142, 1276.54960008, 990.571280057, 196.88970278,          398.206798139, 62.7937319455, 111.92795268), Rdp1D.sense = c(197.403527744,          168.849473212, 399.588620598, 68.0531849874, 128.833494553,          30.8082175235, 59.9086910765, 134.404417978, 24.2425410143,          85.4825519212), Rdp1D.antisense = c(86.097230688, 254.128565899,          388.725581635, 846.769716459, 82.1986385122, 281.872704472,          49.97022677, 77.2892621321, 44.6799202033, 1.60870068737),          Wt.sense = c(150.835381912, 132.061554165, 607.58955888,          65.8027665102, 89.3919476073, 83.4968237124, 7.90112304898,          10.714546021, 5e-04, 5e-04), Wt.antisense = c(150.374084859,          131.8668254, 659.887826114, 65.7197527173, 45.4289405873,          40.4019469576, 7.40733410843, 8.83958796731, 43.5756796108,          12.3289419357), Rdp1D.Rrp6D.sense = c(278.940777843, 227.050371919,          266.352999304, 43.8265653895, 86.2348572529, 5.1007112686,          63.5315969071, 138.590379851, 17.1377883364, 47.2571674648         ), Rdp1D.Rrp6D.antisense = c(122.812370852, 165.478532861,          262.217884557, 315.685821866, 196.899101029, 181.217276367,          64.9492021228, 111.77461648, 62.2771817975, 20.3596716974         ), Dcr1D.sense = c(5e-04, 120.491414743, 1325.93762159, 546.346320658,          5e-04, 5e-04, 66.3486618734, 5e-04, 5e-04, 5e-04), Dcr1D.antisense = c(5e-04,          8346.5035927, 1479.42139464, 37845.8172699, 5e-04, 28845.1503745,          1194.26663745, 5e-04, 647.428121154, 5e-04), Er1D.sense = c(387.657094655,          332.176880363, 570.413411676, 136.333361806, 228.023187499,          5e-04, 24.0778502632, 62.6341480521, 32.1717485621, 5e-04         ), Er1D.antisense = c(382.664804454, 343.714717963, 618.13806355,          205.325286003, 162.81296098, 145.575708252, 15.3360737154,          30.5382985528, 5e-04, 13.8803856753), Rrp6D.sense = c(716.001844534,          605.02996247, 444.912126049, 213.265421331, 398.7252034,          73.8307932225, 90.5802807096, 172.093792998, 5e-04, 135.365316918         ), Rrp6D.antisense = c(690.534019176, 592.944889017, 409.413915909,          247.869927895, 160.655498164, 371.504850116, 56.7600331059,          119.421944835, 16.7787329876, 20.0208426702), Mlo3D.Ago1D.sense = c(119.466474712,          329.741829677, 993.941348153, 1072.99933641, 5e-04, 377.539482989,          113.878508361, 50.428609435, 5e-04, 5e-04), Mlo3D.Ago1D.antisense = c(120.543892198,          2711.8968975, 1257.1652648, 11870.674213, 125.725150183,          8902.64920707, 206.72008398, 37.8215820763, 5e-04, 5e-04),          Ago1D.Clr3D.sense = c(184.712264891, 179.831117561, 444.487152139,          162.69482267, 202.293495599, 5.61159966339, 63.6233691066,          90.544306737, 5e-04, 170.284591079), Ago1D.Clr3D.antisense = c(57.5740294693,          67.5638155026, 386.644572497, 102.906975334, 79.4664091704,          2.1204925561, 14.4184581702, 35.3125846275, 5e-04, 5e-04),          Dcr1D.Rrp6D.sense = c(45.8846113251, 63.7325750806, 360.192351832,          126.841847799, 277.614908589, 54.2822292313, 33.9452752392,          83.1313557186, 5e-04, 12.8242338794), Dcr1D.Rrp6D.antisense = c(19.3160147626,          55.5834301591, 363.594792664, 183.776577157, 18.3768674716,          322.564097746, 17.907465048, 33.1088927537, 5e-04, 5e-04),          Ago1D.sense = c(29.0628360487, 31.9691923002, 387.82120669,          42.2593617334, 64.0004397647, 68.0567121551, 65.0088334947,          189.345502766, 5e-04, 26.5639424914), Ago1D.antisense = c(10.918535798,          84.6095118936, 373.635073395, 345.064708329, 40.1150042497,          266.756186351, 4.38085691952, 5e-04, 5e-04, 5e-04), Mlo3D.sense = c(2798.34040679,          2353.07409522, 330.364494647, 781.101862885, 1312.81871554,          376.811874795, 124.564566466, 353.76677093, 5e-04, 31.5118039429         ), Mlo3D.antisense = c(2532.2553647, 2248.78653802, 292.881120203,          1246.84984213, 1981.14439149, 564.070923014, 164.753382721,          449.669663275, 5e-04, 5e-04), Ago1D.Rrp6D.sense = c(86.379996345,          90.4014346003, 468.105009795, 104.668452639, 203.155350014,          62.3955638527, 44.5603393841, 84.3076975857, 16.0419716595,          42.5345756816), Ago1D.Rrp6D.antisense = c(45.0506816078,          80.7182081997, 481.700138654, 206.646370214, 67.1332741403,          129.669542952, 23.7209335341, 26.0270063646, 28.9823086155,          16.4901597751)), .Names = c("chromosome", "start", "end",      "strand", "annotation", "Abp1D.sense", "Abp1D.antisense", "Iki3D.sense",      "Iki3D.antisense", "Iki3D.Rrp6D.sense", "Iki3D.Rrp6D.antisense",      "Rdp1D.sense", "Rdp1D.antisense", "Wt.sense", "Wt.antisense",      "Rdp1D.Rrp6D.sense", "Rdp1D.Rrp6D.antisense", "Dcr1D.sense",      "Dcr1D.antisense", "Er1D.sense", "Er1D.antisense", "Rrp6D.sense",      "Rrp6D.antisense", "Mlo3D.Ago1D.sense", "Mlo3D.Ago1D.antisense",      "Ago1D.Clr3D.sense", "Ago1D.Clr3D.antisense", "Dcr1D.Rrp6D.sense",      "Dcr1D.Rrp6D.antisense", "Ago1D.sense", "Ago1D.antisense", "Mlo3D.sense",      "Mlo3D.antisense", "Ago1D.Rrp6D.sense", "Ago1D.Rrp6D.antisense"     ), row.names = c(NA, 10L), class = "data.frame")      The question asked above is when you have a `data.frame` with all the data. What if I now want to create a `list` so that each entry in the list is actually a `data.frame` with a structure similar to one given above. How do I combine the boxplots in the lattice? Does the `ggplot2` have a solution for this? Can someone guide me to such a solution?


title: Inaccurate R prediction

text: After development of recommendation engine with the R, before removal of outliers from data-set value of residual standard error was 1351 and after removal of outlier its 656. Still there is no accurate prediction which gives 10% correct(near) prediction. For more fitting i also have tried polynomial model with two ,three and four degree but still no improvement. Is there any most important thing to consider without R-squared or adjusted R-squared. Where i am using dataset with linear regression model for prediction of product purchase revenue on the base of total numbers of time product added to cart, removed from cart, total numbers of page views of product page. For checking model prediction accuracy i am considering only minimum residual standard error. Here is Model summary               > summary(model_out)          Call:     lm(formula = yitemrevenue_out ~ xcartaddtotalrs_out + xcartremove_out +          xproductviews_out + xuniqprodview_out + xprodviewinrs_out,          data = as)          Residuals:         Min      1Q  Median      3Q     Max      -2671.1  -173.6   -83.4   -42.9 14288.6           Coefficients:                           Estimate Std. Error t value Pr(>|t|)         (Intercept)          3.992e+01  1.254e+01   3.183  0.00147 **      xcartaddtotalrs_out -7.888e-03  2.570e-03  -3.070  0.00216 **      xcartremove_out     -3.410e+01  2.431e+01  -1.403  0.16076         xproductviews_out    1.248e+01  1.222e+00  10.215  < 2e-16 ***     xuniqprodview_out   -1.350e+01  1.487e+00  -9.076  < 2e-16 ***     xprodviewinrs_out    3.705e-04  5.151e-05   7.193 7.62e-13 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 656.4 on 3721 degrees of freedom     Multiple R-squared: 0.1398, Adjusted R-squared: 0.1386      F-statistic: 120.9 on 5 and 3721 DF,  p-value: < 2.2e-16       Thanks


title: sampling method in moderation analysis

text: correct me if I'm wrong, I've read that my sample should be normally distributed if I'm going to use moderation analysis. But if in case I'll be using purposive sampling (I'm following certain criteria / conditions )what can I do as a solution in the problem of normal distribution? my last question will be, is moderating variables needed to be tested with normal distribution aside from Independent variable? Thank you very much for the response.


title: First differences interpretation

text: I have a small issue regarding first difference models. Suppose we have $$ \Delta y_{t} = \beta_{0} + \beta_{1} \Delta \log{y_{t-1}} + \Delta u_{t} $$ The interpretation of the delta would be the same as if I did with: $$ y_{t} = \beta_{0} + \beta_{1} \log{y_{t-1}} + u_{t} $$ That is correct right? However if we have: $$ \Delta y_{t} = \beta_{0} + \beta_{1} \log{p} + u_{t} $$ Then the interpretation would be if p increases by 1%, the increase in $y_{t}$ would be $\beta_{1}$ right? Thanks in advance for clearing things up!


title: Find the pdf of Y when pdf of X is given

text: $$ f_{X}(x) = \frac{3}{8}(x+1)^{2} ,\ -1 < x < 1 $$ $$Y = \begin{cases} 1 - X^{2} & X \leq 0,\\ 1- X, & X > 0.\end{cases}$$ I started with : $$ F_{Y}(y) = 1 - P(Y \leq y) $$ $$ = 1 - [P(-(1-y)^\frac {1}{2} < X < (1-y)] $$ From here, I can get $F_{Y}(y)$, and differentiating it will give me $f_{x}(x)$. But the answer I am getting for pdf is not the desired answer. Am I doing anything wrong? Thanks for your help.


title: SPSS three way dummy variables and descriptive statistics

text: I have two dummy variables, X1 and X2. Each one is coded into 0 or 1. Here is the breakdown:                   alcohol use                  X7 X8         none     0   0         moderate 1   0         severe   0   1      I don't understand how to tell SPSS that when X7 is 0, AND X8 is 1, that means 'Severe'. The problem arises when I try to do descriptive statistics. It gives me a mean of the 1's and 0's, which is a useless value of 0.5. How are descriptive statistics shown with dummy variables, and especially when the dummy variable is coded across two separate variables?


title: Interpreting out of bag error estimate for RandomForestRegressor

text: I am using RandomForest regressor on my data and I could see that the oob score was obtained to be 0.83. I am not sure how it came out to be like this. I mean my targets are high values in the range of 10^7. So if it's MSE then it should have been much higher. I don't understand what 0.83 signify here. I am using python's RandomForestRegressor of the sklearn toolkit. I do model = RandomForestRegressor(max_depth=7, n_estimators=100, oob_score=True, n_jobs=-1) model.fit(trainX, trainY ) Then I see model.oob_score_ and I get values like 0.83809026152005295


title: Learning R Programming Question

text: I am trying to learn R programming and I am having trouble with a problem. I am trying to set a probability of success $p$ and find out what the chance of success is for $n$ trials. Basically a Bernoulli trial but I am trying to implement it myself using a `while` loop to learn R better. I know that the loop structure in R is similar to many languages, i.e., `while (condition) {body}` but not sure how to go about accomplishing my task.


title: How to interpret chi-squared test for main effect in presence of interaction term

text: I have run a sequence of nested multi-level models using lmer, and looked at both the summary output giving all the coefficients, etc, as well as the Anova output which does a chi-squared test on each of the factors. In the first model without the interaction, the factor medium is significant, looking at both the chi-squared test output and at the coefficients from the summary output. In the second model when the interaction is added, the first-level factor medium is no longer significant in the chi-squared test output yielded by the Anova function. (It isn't significant in the summary output where we get the coefficient information, but I believe I understand how to interpret that output - it just means that the medium is not significant for the reference group I believe?) But how do I interpret the non-significant outcome of the chi-squared test for the first-level factor medium in the presence of the added interaction? If it helps, here is a simplified version of my model:               success~medium+gender+gender*medium      And here is the Anova (chi-squared tests) output for Model 1 (without the interaction) and Model 2 (with the interaction): Model 1:                               Chisq           Df  Pr(>Chisq)     medium          42.03818323     1   0.0000     gender          6.082823479     1   0.0137      Model 2:                               Chisq           Df  Pr(>Chisq)     medium          2.23141795      1   0.1352     gender          12.95254908     1   0.0003     medium:gender   6.920056392     1   0.0085      I'd be happy to read a reference on this, but I just can't find anything in any of my books on how to interpret the chi-squared tests in the presence of an interaction - they all seem to focus entirely on interpreting the regression coefficients only. So if anyone could help to explain the correct interpretation, or refer me to a good reference that specifically talks about how to interpret the chi-squared tests of main effects in the presence of interactions, I'd be grateful!


title: Finding "weak ties" in network data

text: So I am working on a new project looking at formal and informal networks between businesses in the same industry. Namely, I am looking at joint ventures, fractional acquisitions, minority share ownership, and board of director relationships between companies. I want to see if I can identify the occurrences of "weak ties" (Granovetter, 1973), meaning ties between densely interconnected clusters. However, I was not sure if software was available to do that type of cluster and then weak tie identification? I am pretty familiar with the ERGM models and the Statnet packages, but I have not seen an functionality for this in those packages. I might be able to use a latent space model (Handcock, 2006) to identify clusters and then eyeball the weak ties, but that seemed a bit arbitrary--though I might have to resort to that. Does anyone know of any better software packages or code for cluster/community detection and then the identification of weak ties?


title: Why aren't all tests scored via item analysis/response theory?

text: Is there a statistical reason why item analysis/response theory isn't more widely applied? For instance, if a teacher gives a 25 question multiple choice test and finds that 10 questions were answered correctly by everyone, 10 questions were answered by a really low fraction (say 10%) and the remaining 5 were answered by roughly 50% of people. Doesn't it make sense to reweight the scores so that hard questions are given more weight? And yet, in the real world tests almost always have all questions weighted equally. Why? The below link discusses discrimination indices and other measures of difficulties for choosing which questions are best: http://fcit.usf.edu/assessment/selected/responsec.html It seems though that the method of figuring out the discrimination index of questions is only used in a forward-looking way (eg., if a question doesn't discriminate well, toss it). Why aren't tests re-weighted for the current population?


title: Wald $\chi^2$ value for a probit model

text: I have fit a probit model in R. However, in addition to the $z$- and $p$-values, I would like to know the Wald $\chi^2$ values for the individual explanatory variables. How do I obtain / calculate the Wald $\chi^2$ values for parameter estimates?


title: Estimation the standard error of correlated (binomial) variables

text: Say we have N voters from some group (N binomial variables). Each voter has a preference on voting about some issue (probability of "success" p), which depends on the preferences of other voters. Given a sample of votes for each voter (an estimation of their preference on that issue), and an estimation of the correlation matrix between each two voters. Can we give a better estimation of the standard error of the average of the proportions (than the one offered by the standard $\sqrt{\frac{p(1-p)}{n}}$?


title: Does the principle of indifference apply to the Borel-Kolmogorov paradox?

text: Consider Jaynes' solution to the Bertrand paradox using the principle of indifference. Why doesn't a similar argument apply to the Borel-Kolmogorov paradox? Is there something wrong with arguing that since the problem does not specify an orientation for the sphere, rotating the sphere should not affect the resulting distribution arrived at by the chosen limiting process?


title: Proportional hazard ratio formula in R

text: I'm new proportional hazard ratio models and new-ish to R, so I suspect this is a basic question: The scenario is modeling the speed of an effect through hazard ratios. People are joining an organization, and most are invited by someone else within the organization*. The speed measured is the time between when the invitation is made and the person joins the organization, called the `registration_interval`. The way that seems evident to model this in R is like so:               analysis<-coxph(Surv(registration_interval) ~ Factor1 + Factor2 + Factor3 ...)      However, there's also the prospect that people are faster/slower to join the organization as the date of their invitation moves along (ex. people invited in 2012 might join faster than people joining in 2009), and this could interact with the various factors. What is a sensible way to incorporate this `date` information? *We only have data for those that do ultimately join the organization. We consider our model as simply modeling "joiners", and not right-censoring everyone whose invitation-registration intervals are infinitely long. Edit for clarification: We know when people are invited, so on any given date there are only a certain number of people invited who could join. If everyone were invited all at once the proper model setup would just be:               analysis<-coxph(Surv(registration_interval) ~ date + Factor1 + Factor2 + Factor3 ...)      But since the `registration_interval` is determined by two dates, there is some interaction between that and `date`. The hazard function would, for example, look quite different if there is a span of time in which nobody is invited, versus a span of time in which thousands of people are invited. Perhaps this is all packaged within `Surv`, and I'll look into it. Edit regarding putting both `date` and `registration_interval` or join time information in Surv: Nope. That appears to be a red herring.


title: Build a covariance matrix for GAM

text: I am using a GAM in R to compare time series data from 3 countries. The data sets are of hourly measurements for one year. The main aim here is to show when at which time of the day and day of the year are the data well matched to the mean series. The script below shows my attempt:               ## Ozone measurements for three countries in Europe     ## Find similarities between time series     require(plyr)     require(lattice)     require(mgcv)          TopFolder <- list("http://www.nilu.no/projects/ccc/onlinedata/ozone/CZ03_2009.dat"     ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/CY02_2009.dat"     ,"http://www.nilu.no/projects/ccc/onlinedata/ozone/BE35_2009.dat"     )          ## create variable for the data      data = ldply(TopFolder, header = TRUE, read.table, sep = "", skip = 3)          ## define Ozone levels     Ozone <- data$Value     Ozone[Ozone==-999] <- NA     Ozone <- data.frame(Ozone)          ## define Datetime -  need to concatenate arrays     DateTime <- paste(data$Date,data$Hour, sep = " ")     Date <- as.POSIXct(DateTime, format = "%d.%m.%Y %H:%M")          ## define countries      Countries <- c("Czech","Cyprus","Belgium")     Country <- data.frame(Country = rep(Countries, each = 8760))          ## bind together     Dat <- cbind(Ozone, Country = Country)     Dat <- transform(Dat, Doy = as.numeric(format(Date,format = "%j")),                          Tod = as.numeric(format(Date,format = "%H")),                    DecTime = rep(seq(1,365, length = 8760),by = 3))          ## plot the Ozone data     xyplot(Ozone~DecTime | Country, data = Dat, type  = "l", col = 1,            strip = function(bg = 'white',...)strip.default(bg = 'white',...))          ## generalised additive model     mod1 <- gam(Ozone ~ Country + s(Doy,bs = "cc", k = 20) + s(Doy, by = Country, bs = "cc", k = 20) +        s(Tod,bs = "cc", k=7) + s(Tod,by = Country,bs = "cr", k=7),                  data = Dat, method = "ML")          plot(mod1,pages=1,scale=0,shade=TRUE)          xyplot(resid(mod1) ~ Doy | Country, data = Dat, type = c("l","smooth"))          mod2 <- gamm(Ozone ~ Country + s(Doy,bs = "cc", k = 20) + s(Doy, by = Country, bs = "cc", k = 20) +        s(Tod,bs = "cc", k=7) + s(Tod,by = Country,bs = "cr", k=7),                  data = Dat, method = "ML",correlation = corAR1(form = ~ DecTime | Country))      One problem of this is that I am using an additive model with correlated errors i.e. in each time series of hourly measurements, a high concentration is followed by another high concentration. This can be addressed by adding a correlation matrix into the GAM: However, R throws out an error as it is unable to get any more RAM from the OS. So, to tackle this I need to: (1) from mod1 - fit a time series model (mod3) to the residuals (2) obtain the acf(1) from mod3 (2) use the acf(1) to build a covariance matrix that we can then pass into mod2 thus simplifying the correlation matrix. Any advice on how to to complete some or all of the three steps above would be highly appreciated. Although I think the first point that I an unsure of is how to extract the residuals for the different countries in mod1?


title: Dimension reduction techniques for very small sample sizes

text: I have 21 socio-economic and attitudinal macro-level variables (such as percentage of mothers aged 24-54 not employed, percentage of children aged 3-5 in nursery schools and so on). I also have data on the proportions of grandparents which provided intensive childcare. Most of the socio-economic variables which I selected are highly correlated with childcare provision (for instance, there is a negative correlation between proportion of mothers employed part-time and provision of grandparental childcare). Ideally, I would like to create a typology of different kinds of countries. My hope would be to use some kind of dimension reduction technique whose components or factors would make some intuitive sense (e.g. attitudes towards family and gender, labour market structure, family policies). Or, alternatively, assess which of the 21 macro-level indicators best explain the variability in childcare provision across countries. My main problem is that I only have 12 European countries. I reckon that PCA and factor analyses are not appropriate techniques with so few cases. Am I correct? I was told to try use qualitative comparative analysis or multiple correspondence analysis, although to my understanding the latter techniques are more appropriate for binary (or categorical) macro-level indicators (whereas mine are percentages or continuous variables).


title: On cophenetic correlation for dendrogram clustering

text: Consider the context of a dendrogram clustering. Let us call _original dissimilarities_ the distances between the individuals. After constructing the dendrogram we define the _cophenetic dissimilarity_ between two individuals as the distance between the clusters to which these individuals belong. Some people consider that the correlation between the original dissimilarities and the cophenetic dissimilarities (called _cophenetic correlation_ ) is a "suitability index" of the classification. This sounds totally puzzling to me. My objection does not rely on the particular choice of the Pearson correlation, but on the general idea that any link between the original dissimilarities and the cophenetic dissimilarities could be related to the suitability of the classification. Do you agree with me, or could you present some argument supporting the use of the cophenetic correlation as a suitability index for the dendrogram classification ?


title: Classification with GBM in R and imbalanced class sizes

text: I'm dealing with a supervised binary classification issue. I'd like to use the GBM package to classify individuals as uninfected/infected. I have 15 times more uninfected than infected individuals. I was wondering if GBM models suffer in the case of imbalanced class sizes? I didn't find any references answering this question. I tried to adjust the weights by assigning a weight of 1 to the uninfected individuals and a weight of 15 to the infected, but I obtained poor results.


title: SVM kernel mapping, finding boundaries in projected space

text: I have a question about the support vector machine (SVM) kernel trick. How do you find the boundaries of the training data set in kernel projected space? Is that the same boundaries as you can obtain in original input space? Does kernel mapping change the data distribution? I am new to SVM, can someone please give me a straightforward explanation or links to related information?


title: Latent groups in Bayes net with BUGS

text: I'm modeling a Bayes net with OpenBUGS, and I find problems to specify some of the parameters and their priors. The aim of the model is to identify latent groups in the data from a sample of human participants. Eventually, I need a belief distribution over the range of 5 possible latent groups, for each individual, conditional on the observed data. In my model, the observed data `D` come from a normal distribution with mean `theta`. The value of `theta` is fixed and depends on the state of a third variable, `G`, which represents the latent groups. This works such that, if the variable `G` possesses value 1, `theta` equals `x`, if `G=2`, `theta=y`, and so on, up to an arbitrary and low number of predefined latent groups (it's actually 5 groups in my current problem). Note that the number of latent groups and their respective values for theta are defined by the experimenter, so this is not a typical clustering problem. In fact, each latent group corresponds to the prediction of one different theory about the experimental situation. I don't know how to include the specification for `theta` in the BUGS model. That should be added in the code below (I simplified it a lot for posting here):                   model {         for (i in 1:NumberOfParticipants){             D[i] ~ dnorm(theta[i],lambda[i])            lambda[i] ~ dgamma(.01,.01)           p<-1/5           G[i] ~ dcat(p[]) # assume that all latent groups are equiprobable a priori?           }          }       I need to indicate in this code how a given value of `theta` is determined by the latent group to which the current participant silently belongs. I admit that my approach to the problem could be completely wrong.


title: Probability distribution for varying probabilities in R

text: I'd like to plot the probability distribution for a set of binomial trials in R, the catch is that each trial has an independent probability of success (which I have in vector form). So, in R if I were doing this with coin tosses: `plot(x,dbinom(x,10,.5))` would work, where `x` is `0:10`. This shows the probability distribution by plotting number of successful trials on the x-axis with percent of the time that specific results is achieved as the `y` (4 successes is 20.5%). That said, how do I plot the same graph for 10 discrete trials with varying probabilities. For instance if `odds<-c(.1,.8,.2,.2,.3,.7,.9,.99,.05,.5)`, was the probability of success for each independent toss, and I wanted to see a probability distribution?


title: Is there a distribution appropriate for a continuous variable skewed toward zero and able to include zero?

text: I am interested in modelling the impact of some environmental parameters on a concentration of measured phytoplankton pigment. The concentration of pigment is skewed so that low concentrations are more frequent than high ones and there are occasions of 0 (zero) concentration. I was hoping to use a GLM model with a Gamma distribution, but the zero values prevent such a fitting. My question is if there are other options for me. I thought to simply add a very small number to all concentrations (e.g. 0.00001) and proceed, but perhaps this is ill-advised. Given that the response variable is continuous, I am not sure what my options are - Poisson and Negative Binomial distributions are only for count data? Any advise would be greatly appreciated.


title: "Wrong Sign" On Regression Coefficients - Hierarchical Multiple Linear Regression

text: I am analyzing my data on the relationship between spirituality and negative emotional states (depression, anxiety, and stress) using a hierarchical multiple linear regression. Everything seemed to be "okay" with the results until I checked the signs of the regression coefficients!! The results are showing a positive relationship between depression and a spirituality factor which includes having a reason for living. Virtually all of the research in this area shows that there should be a negative relationship between these two factors, i.e. that higher levels of depression are correlated with lower levels of having a reason for living (without any interaction terms or moderating variables). Here are the steps that I have completed so far in this three-part data analysis (1st data analysis - Principal components analysis / 2nd Part - EFA (Exploratory Factor Analysis) / 2nd data analysis - hierarchical multiple linear regression):   1. Created spirituality questionnaire and collected data from sample (n = 189).   2. Used a principal components analysis to identify how many factors to extract in EFA.   3. Derived four factors from EFA.   4. Factors included log transformed questions which changed the directionality of the questionnaire (i.e. questions which were "positively" worded became negatively worded).   5. All questions on questionnaire were standardized (transformed to z-scores) to correct for directionality problem.   6. Negative emotional states were examined for normality. Two of the three negative emotional states had to be log transformed.   7. Ran 3 hierarchical multiple linear regressions with DV = Either log depression, log anxiety, or stress and IVs = Age (first "block"), four spirituality factors (second "block") using enter method.   8. Noticed that signs were incorrect, so I compared log transformed versions of DVs to non-transformed versions AND standardized DV (z-score) versions (see this article: The Independent Sign Bias: Gaining Insight from Multiple Linear Regression [PDF] --> pg. 1)   9. IV having to do with meaning in life had correct sign in non-transformed version of depression HOWEVER it looks like there is a collinearity issue between two of the spirituality IVs with the non-transformed version which may affect the signs!!! (Condition index > 15 and Variance proportions: .81 and .40)   10. Looked at bivariate correlations to see if the signs were correct, and oddly enough the signs were accurate (this was only true for the non-transformed versions of the hierarchical multiple linear regressions). None of the correlations were above .6.   11. Removed Age as an IV to see if that made a difference. It did not.    12. Removed spirituality IVs one at a time with stepwise method rather than enter method and this only increased the amount of collinearity between the variables that remained in the final model. I would greatly appreciate any suggestions on how to correct for the wrong signs and/or if I did anything incorrectly in my three-part data analysis.


title: Classification model that allows changing priors at prediction time

text: I would like to know if there's a way to build a classification model in R that would allow me to change the class weights at prediction time. The scenario where I would want to do this: I have a group of products, that are one of three different types (A, B, and C). The types occur in in observed frequencies of 50% A, 30% B, 20% C. I train a classification model to recognize the type based on some set of attributes common to all three types. The classification model is not perfect, but has decent accuracy. If I have some new observations of a product, I can run them through the classification model and get a predicted type (or actually probabilties of being any of type A, B, or C). After the model is trained, I discover that each product was produced from one of 10 different machines (M1 ... M10). Each machine is capable of producing more than one type of item, but does so with a specific set of probabilities. e.g.   * M1 = (50% A, 25% B, 25% C)   * M2 = (33% A, 33% B, 34% C)   * M3 = ( 0% A, 50% B, 50% C)   * M4 = (10% A, 85% B, 5% C) The machines produce items at different rates that contribute to stable overall population frequencies of 50/30/20. If I don't know which machine produced a newly observed product, I can continue using the existing model as before. But if I also know that the new product was produced by, say, machine M4, I would like to be able to adjust the model's posterior probabilities to account for M1's class priors. Conceptually, something like:               predict(model, newdata=new.observation, new.class.weights=c(0.1, 0.85, 0.05))      Is there an R package that would allow me to build a classification model where I can specify a different set of priors at the time of prediction, rather than when the model is built? Or is there a way to adjust the probabilities from the model to account for the newly known prior probabilities? I'm open to any classification algorithm or R package, if it matters.


title: Estimate multivariate distribution of dependent variables

text: Let $H$~$\mathcal{N}(\mu_H,{\sigma_H}^2)$ and $S$~$\mathcal{N}(\mu_S,{\sigma_S}^2)$. $H$ and $S$ are dependent of each other and their covariances are unknown. Let $G = H - S$ and the mean and variance of $G$~$\mathcal{N}(\mu_G,{\sigma_G}^2)$ are known. Is there a way to I estimate the multivariate normal distribution of the vector $(H,S)$? Sorry to bother, if this question is trivial, but I'd appreciate help. JBashir


title: Power of t-test for aggregate data

text: I have data that have been aggregated with unknown variance into a number of composite averages. These data are generated from pooled assays, such as HIV viral load in a composite of 10 vials of blood, or cystic fibrosis protein concentration in a composite of 10 0.01cc biopsies of lung tissue. We have treated and controlled groups and are interested in testing whether there is a difference. We must calculate power based on preliminary results (known effect size/hypothesized mean difference) to justify assay costs and find the most cost-effective aggregation level (e.g. 10 samples per assay vs. 20 samples per assay). For the primary analysis, we will just use a t-test to determine if there are differences in the endpoint between groups. As I said, we do not have an estimate of the individual level variance in each composite, so inverse variance weighting is not an option. However, we do know how many samples were in the composite, so we can use frequency weighting. How do I calculate the non-centrality parameter for the t-test test statistic with unknown variance? Also, what are the effective degrees of freedom? Note simulation is not an option, I need the exact distribution of the test statistic in order to use optimization routines to find maximum power under a fixed cost.


title: How to estimate the probability that the mean of an unknown distribution is over a threshold given small sample size

text: I am trying to quantify my concerns regarding a proposed incinerator in our community. The company is basing its potential to emit dioxins (a class of chlorinated organic compounds with a reference dose of 1.7E-8 grams per 150 lb person per year) on three emission measurements (mg/m^3 exhaust): 0.0002139; 0.0000014; 0.00000186 The company has submitted their potential to emit as the mean of these three values, which is near legal limits. If the plant is built, they may be required to measure dioxins only once on completion, or possibly every year thereafter. My initial thoughts have been to find the probability that these samples were drawn from a distribution with a mean above a certain threshold (e.g. the legal limit). As sample measurements, I would assume a normal distribution, though I guess with some non-zero skewness. I don't have any prior knowledge about the variance other than the samples. But, I am now confronted with an infinite number of possible distributions, and am unsure how to proceed. I suspect this is a well studied area of statistics (risk analysis?), and any pointers in common methods would be much appreciated. I am a biochemist, and statistics is not my first language.


title: Can you determine if variances are equal when quartiles are similar?

text: I am comparing 2 different data sets. I have broken down each set into quartiles. I am looking at the percentage of increase between the quartiles (example there was a 125 percent increase between the 2nd and 3rd quartile for one data set and a 123 percent increase between the 2nd and 3rd increase of the other data set) can I say that the variability is the same for both groups if the percentage of increase between all the quartiles are about the same between the two datasets?


title: How should I compare two methods for obtaining the frequency of several events?

text: I am using two different methods to obtain the frequency of several events. Can I use a scatter plot to graph the frequency of an event obtained by Method A against the one obtained by Method B? And afterwards, report the correlation coefficient to state that, since they strongly correlate, they can be used interchangeably depending on the resources available? Or should I use a Bland-Altman plot or a least products regression analysis?


title: R codes for variation of information criterion using "mclust"

text: I am developing model-based clustering. First, I developed model-based clustering in R using "mclust." Next, I wanted to take 75% of the sample, re-run model-based clustering and compare the results with the results from the entire dataset using variation of information or rand index. However, I am getting stuck with the codes. Here are the codes made available for variation of information on CRAN.               cl1 <-sample(1:30, 10, replace=TRUE)     cl2 <- c(cl1[1:5], sample(1:3, 5, replace=TRUE))     vi.dist(cl1,cl2)     vi.dist(cl1,cl2, parts=TRUE)      Here is how I developed my clustering solutions:               model <-Mclust(data[,18:22])      model based clustering solution using entire dataset and data is the name of my dataset using columns 18 to 22               random <- data[rbinom(nrow(data), 1,.75)==1,]      developed a random sample               randommodel <-Mclust(random[,18:22])      model-based clustering using the random sample               cl1 <- (model$classification)     cl2 <- (randommodel$classification)     vi.dist(cl1,cl2)      my attempt at variation of information using R codes, result failed because `cl1 and cl2 are not the same length`. So, how do I make the two solutions the same length given that the two solutions have a different number of observations? Or, am I trying to use variation of information incorrectly?


title: Workflow for using R with LyX for a statistical analysis?

text: For some time, I wanted to stop copy-pasting my R results into word, but climbing the LaTex mountain seemed to much to be worth it. Recently, I came to discover LyX, as a laymen's solution to people like me who do not wish to code their text, but do wish to combine R analysis with text. However, I found there is very little, **updated** , documentation about LyX+R - which leads me to writing my question:   1. What tools do you combine in your R+LyX workflow? (do you combine text editors on top of LyX? which? and why?)   2. What is your common folder structure for an analysis project?   3. What is the order of steps you take for constructing your analysis? What code do you keep in file.r? and what in other types of files? (images, .RData backups, .TeX, and so on)   4. Do you use different work strategies for different types of project (due to size of project or size of dataset)   5. What R packages do you combine in your work? and how?   6. Recommended links? Thank you for any input!


title: Measuring association using Cramér's $V$

text: I had measured the association between two categorical variables some time back. Now, I have even more data compared to the last time. I find that the Cramér's _V_ values has decreased now. Is this because of the increase in the number of data points? Can we at all attribute this change in Cramér's _V_ because of the increase in data points? What can be the reason behind increase/decrease in Cramér's _V_ value?


title: Does and 95% CI containing 1 imply no difference in IRR calculations in the same was as in OR calculations?

text: My understanding is that with an odds ratio (OR) that a 95% CI containing 1 means that there is no difference between the odds of something occurring between groups. I wondered if this holds true for incident rate ratios (IRR)? Many thanks, Ewen


title: Kappa statistic in plain English?

text: I'm reading a data mining book and it mentioned the Kappa statistic as a mean for evaluating the predication performance of classifiers. However I just can't understand this. I also checked wikipedia but it didn't help too. How does Kappa help in evaluating the predication performance of classifiers? What does it tell?!! I understand that 100% kappa means the classifier is in total agreement with a random classifier, but I don't understand how does this help in evaluating the performance of the classifier! What does 40% kappa mean? Does it mean that 40% of the time, the classifier is in agreement with the random classifier? If so, what does that tell me or help me in evaluating the classifier?!!


title: Libraries for (label propagation algorithms/frequent subgraph mining) for graphs in R

text: **General description of the problem** There is a graph where some of the nodes have a certain type(There are about 3-4 types). For other nodes, type is not known. I want to predict, based on my graph, for the nodes with unknown type their "most probable" type. **Possible framework** I guess the general framework for such tasks is called `label propagation algorithm`. according to literature on the topic. Here are some examples: one, two Another often mentioned topic is `Frequent Subgraph Mining`, which includes algorithms like `SUBDUE`,`SLEUTH`, and `gSpan`. **Found in R** The only label propagation algorithm I managed to find in `R` is `label.propagation.community()` from `igraph` library. However as the name suggests it is mostly for finding communities, not classification of nodes. There also seems to be several references to `subgraphMining` library, (here for example) but looks like it is missing from CRAN. **Question** If someone could suggest libraries/frameworks for the task described, I would be very grateful. Made the question as specific as I could, hope that shall do)


title: Good papers with reproducible analysis requiring only the basics

text: I'm looking for papers or other examples of research where the statistical analysis done would be within the grasp of someone who has done an introductory stats course. Ideally the datasets would also be available online. The idea is to pass these to students and examples they can follow and reproduce.


title: Reporting Fixed Effects as (partial) correlations?

text: I'm doing a linear mixed effects analysis in which I'm really only interested in one of the fixed effects. I have several other fixed effects and a random intercept term, but none of them are relevant, nor is the overall goodness-of- fit of the model (other than the fact that adding my regressor of interest (ROI) improves it significantly). So, I now want to report the size of my interesting fixed effect in an intuitive way. In a simple multiple regression analysis, I would compute the partial correlation between my ROI and the outcome measure by regressing out all other factors from both the ROI and the outcome measure, and then computing the correlation between the residuals from these two "nuisance regressions". Or, equivalently, taking the beta coefficient on the ROI and multiplying it by the ratio of the standard deviations of the residuals from the nuisance regressions, i.e. SD(ROI_res)/SD(Outcome_res). I would like to do something similar here (with linear mixed effects), but I have trouble finding the right approach. Certainly, just "copy-pasting" the approach from multiple regression doesn't seem to yield sensible results (as the 2 ways of doing this which were previously equivalent now give different answers). So all of this boils down to 2 questions: 1) Is there a way to obtain partial correlation coefficients for fixed effects in a linear mixed model? 2) If not, then what would be appropriate way to report my effect of interest? Reporting the estimated fixed effect coefficient seems valid but not very insightful in terms of effect size/strength of association (as it is entirely dependent on the scale of the two variables), so my guess is there should be a better way...


title: Missing data (multiple imputation)

text: Can I multiply impute missing data in some of the variables which were for some reason not measured at one time? For instance, for production I have time series data for 1961-2011: however, some of   the predictor variables (rainfall, temperature, humidity and few more) which I suspect could have impact on rice production of Nepal were measured only after 1975. So the data for 1961-1974 are missing! Could multiple imputation help in this case, or is there any other option?


title: Aggregating results from linear model runs R

text: Since regression modeling is often more "art" than science, I often find myself testing many iterations of a regression structure. What are some efficient ways to summarize the information from these multiple model runs in an attempt to find the "best" model? One approach I've used is to put all the models into a list and run `summary()` across that list, but I imagine there are more efficient ways to compare? Sample code & models:               ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)     trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)     group <- gl(2,10,20, labels=c("Ctl","Trt"))     weight <- c(ctl, trt)          lm1 <- lm(weight ~ group)     lm2 <- lm(weight ~ group - 1)     lm3 <- lm(log(weight) ~ group - 1)          #Draw comparisions between models 1 - 3?          models <- list(lm1, lm2, lm3)          lapply(models, summary)


title: Increasing function defined as the covariance bivariate normal

text: Suppose $c>0,\sigma>0$ and $\tau>0$ are fixed real constants. Then I'd like to prove that the function $g_c:(-1,1)\mapsto\mathbb{R}$ defined by \begin{equation} g_c(\rho)=\int_{-\infty}^\infty\int_{-\infty}^\infty\frac{ \{(-c)\vee x \wedge c\} \{(-c)\vee y \wedge c\}}{\sqrt{2\pi(1-\rho^2)}\sigma\tau} e^{\frac{\sigma^2x-2\rho\sigma\tau xy+\tau^2y^2}{2(1-\rho^2)\sigma^2\tau^2}}dxdy \end{equation} is strictly increasing. I have tried to prove by differentiating $g_c$ w.r.t. $\rho$. But, this doesn't help because the result is very ugly. Can anyone help me? Thanks in advance PS: The function $g_c$ can be rewritten as $g_c(\rho)=E( (-c)\vee [(X_{\sigma,\rho}Y_{\tau,\rho})\wedge c])$ for some random variable $(X_{\sigma,\rho},Y_{\tau,\rho})^T\sim N_2(0,\Sigma)$ where $$ \Sigma= \begin{pmatrix} \sigma^2 & \rho\sigma\tau \\ \rho\sigma\tau & \tau^2 \end{pmatrix}. $$ But, I also don't know how to see my problem using this fact.


title: 3d Kolmogorov-Smirnov test

text: Can the KS3D2 test as suggested by Fasano and Franceschini (1987) be used when one of the three variables take discrete values between 0-40? The other two variables are continuous.


title: X and X/20 same distribution?

text: Is it true that if a random variable $$X$$ follows a distribution, say one of these   * Binomial   * Poisson   * Exponential   * Normal   * Uniform Then $$\frac{X}{20}$$ follows the very same distribution as $X$?


title: X is binomially distributed, and Y is binomially distributed. What distribution does X+Y follow?

text: We know that $X$ follows $\mathrm{Bin}(n_1,p_1)$, $Y$ follows $\mathrm{Bin}(n_2,p_2)$, $X$ and $Y$ are independent. What does $X+Y$ follow? I know the answer, that if $p_1=p_2=p$ then $X+Y$ follow $\mathrm{Bin}(n_1+n_2,p)$ but I do not quite get it. I tried to match $X+Y$ to the criterion of Bin distribution:   * Success/fail   * Independence, yes since Y and X are independent, and all $X$'s are independent with themselves, and all Y's are independent with themselves.   * same p. yes, If we are allowed to assume that. Also, what more can we say about $X+Y$?


title: Inter-experimental variation in survival experiment - how to estimate variability?

text: In my animal experiments, I do survival studies, which generate Kaplan-Meier survival curves for each group, which I then compare with an appropriate log rank test. My question is: If I have run a survival experiment with identical variables, say, five times, and the final outcome (in very layman's terms) happens to be slightly different each time, is there a test that I can do estimate the variability (variance) of my experimental runs? One major confounding factor is the fact that survival is a continuous variable over time. Therefore, it is difficult to reduce it to a single statistic (that can then be compared with a test). A lot of people use the median survival (expressed in units of time) as a single statistic surrogate for survival, but it can often be misleading - depending upon the slope of the survival curve, and may not represent the true nature of the survival outcome. Can someone here help? Please also let me know if further clarifications are needed.


title: What technical terms apply to non-linear transforms that retain rank order?

text: I am almost certain there is/are a technical term(s) to describe transformations (e.g $x^2$ or $log(x)$ for positive numbers) which are non- linear transformations of the variable $x$, but retain the rank order of x. But I am drawing a blank and my google searches have be ineffective. What technical terms apply to these sorts of transformations?


title: Mode, Class and Type of R objects

text: I was wondering what are the differences between Mode, Class and Type of R objects? Type of a R object can be obtained by typeof() function, mode by mode(), and class by class(). Also any other similar functions and concepts that I missed? Thanks and regards!


title: Why are statistical properties of mode and median difficult to determine?

text: I have read a book saying that statistical properties of the mode and the median are difficult to determine compared to the mean. I am not entirely sure why is so. Is it correct that because the mean can be put into mathematical functions more easily compared to the mode and median? What does it mean by statistical properties? What does it mean by determining statistical properties? If anyone could provide more explanations or examples of some sort. Thanks for all the help!


title: Box Cox transformation in R

text: I understand how to use the box cox transformation in R and how to get the graph and lambda. These are the things that are confusing me. For simplicity assume this example:                Weight = Gender + Height + Age + Income      gender = categorical variable 1 = male, 0 = female continuous variables - weight, height, age, income. I've done:               model = lm(weight ~ gender + height + age + income)      and applied box-cox to this model which is about 1 so no transformation is needed for weight. The questions are:   1. How do I apply Box-Cox to the 'x' variables to see if they need to be transformed (I read they can be applied to all x variables, the only disadvantage is its time consuming, but this isn't an issue for me).    2. How do I know if they need a transformation for sure e.g. if lambda is 0.4 should I use a square root transformation or does it 'have' to be 0.5? what is lambda is 2?    3. If one or two or all variables need a transformation then how do I adjust the main model formula?


title: Testing "weak" instruments - F test

text: I would like to address the endogeneity problem in my model. Let me go step by step.   I have panel data for 19 countries, 1995-2010. My regression model: `reg gge et (+control variables)` My independent variable is: energy tax (code: et). My dependent variable is: CO2 (equivalents) emissions from industrial processes (code gge). Theory suggests instrumental variables: public debt, public deficit, transport tax Following (Wooldridge 2003) I used:   first differences for my independent and dependent variables in my equation. Than I found instrumental variables: public debt (code: pdeb), public deficit (code: pdef, transport tax (tt). Firstly I would like to test if I have weak or strong instruments. On e.g. slides from Charlie Gibbons (p.8), it says I should check for the inclusion restriction by doing an F test to ensure that this is the case. How should I check for the inclusion restriction to acquire an F test in terms of Stata and SPSS command - simply by putting my IV instead of my dependent variable in a regression or?   e.g.   regression model: `reg gge et` (+control variables)?


title: Shape of posterior distribution

text: I have recently been reading and trying to understand the Bayesian paradigm and looking at various methods that people have been using to estimate the posterior distribution. Now, it seems that most of the intractability comes from the fact that the likelihood terms have a complex form or there are complex integrations to solve. I was wondering what happens when the terms are from some well-defined distribution families. One problem I have been looking as an example has the following setup: So, we want to estimate the parameters for the distribution of a random variable $w$ and we observe some data $d$ and we have a template $x$, so we want to estimate $p(w|d)$. Now assume that the data is observed IID and we have a Gaussian likelihood, so we can write the likelihood as: $$ p(d|x,w,\phi) = \prod_i{\sqrt{\frac{\phi}{2\pi}}\exp^{-0.5(y_i-x_i)\phi(y_i- x_i)}} $$ where $\phi$ is some global noise precision. Also, assume that the prior on $w$ is a multivariate zero mean normal prior $\it{N}(w; 0, \Sigma)$. So, now the posterior distribution is as follows: $$ p(w|d,x,\phi) = \frac{p(d|x,w,\phi)p(w)}{\int p(d)} $$ Now, my question is that since we have a Gaussian likelihood and Gaussian prior, is this posterior tractable and can we say something about how many modes it will have. Also, can someone help me understand intuitively why that may or may not be the case. I would really appreciate any thoughts on this matter.


title: (log) Likelihood at true parameters

text: I am considering the linear regression model \begin{equation} Y_{i} = A t_{i} + B + \varepsilon_{i} \tag{$1$} \end{equation} where $A$ and $B$ are unknown real parameters and $(\varepsilon_{i})_{1 \leq i \leq n}$ are observed values of independent and identically distributed gaussian random variables with mean $0$ and variance $\sigma^{2}$ ($\sigma^{2}$ is unknown). In matrix form, $(1)$ writes : \begin{equation} Y=X\alpha + \varepsilon \tag{2} \end{equation} where $X = \begin{bmatrix} t_{1} & 1 \\ \vdots & \vdots \\ t_{n} & 1 \end{bmatrix}$, $\alpha = \begin{bmatrix} A \\ B \end{bmatrix}$ and $\varepsilon \sim \mathcal{N}(0,\sigma^{2}\mathrm{I}_{n})$. Using Matlab, I generated some data artificially : I took $n=100$, $A=3$, $B=-1$ and $\sigma^{2} = 0.64$ and the parameters estimates I obtain are : $\widehat{A}=3.0096$, $\widehat{B}=-0.9260$ and $\widehat{\sigma^{2}}=0.6337$. The log-likelihood of $(2)$ is the function $$ \ell(A,B,\sigma^{2}) = -n \log(\sigma) - \frac{1}{2\sigma^{2}} \Vert Y-X\alpha \Vert^{2} $$ This may be a silly question but I do not understand why the _true_ values of the parameters do not maximize the log-liklihood in my case. Is it because the estimates I obtained using Matlab are observed values of random variables which asymptotically converge to the true values of the parameters and here, $n$ is not "big enough" ?


title: What econometric model to forecast a seasonal commodity demand while incorporating exogenous Information?

text: I have a monthly commodity demand and try to forecast this series for the next 5 years. Here is a plot: ![Plot](http://i.imgur.com/ccoci40.jpg) Of course, the natural approach to forecast this seasonality would be some kind exponential smoothing (e.g. Holt Winters). But my client requires me to incorporate also information about the state of the economy (e.g. industrial production, household income) as well as other variables (e.g. weather conditions) into the forecast. What I need is some kind of hybrid model, which combines time series forecasting with exogenous information. I thought ARMAX might be an appropiate method (I am working with R) but perhaps there is simpler model. Thank you for any suggestions! edit: Here is the data set in R:               structure(list(Demand = c(7381.08, 6902.08, 7107.08, 5797.08,      4778.08, 4330.08, 4065.08, 3743.08, 4378.08, 4956.08, 5726.08,      6167.08, 8497, 7299, 6565, 4746, 4224, 3755, 3398, 3656, 3634,      4948, 7052, 6268, 7605.75, 8072.75, 6411.75, 5494.75, 3707.75,      3711.75, 3643.75, 3406.75, 4387.75, 5133.75, 6210.75, 6963.75,      9158, 7362, 7772, 5136, 5145, 4394, 3647, 3713, 3965, 5451, 6674,      7419, 7198.42, 7387.42, 7402.42, 5556.42, 4204.42, 4092.42, 3882.42,      3886.42, 4413.42, 5330.42, 7458.42, 7653.42, 7999, 7012, 6063,      5811, 4993, 4299, 3795, 3888, 4395, 5806, 7547, 8249, 8311.42,      6626.42, 6702.42, 5713.42, 4137.42, 3878.42, 3624.42, 3699.42,      5033.42, 6061.42, 7363.42, 8571.42, 8789.08, 8831.08, 6693.08,      6040.08, 5353.08, 4047.08, 3038.08, 3134.08, 3719.08, 6310.08,      8257.08, 9153.08, 9292, 7998, 7566, 6008, 4095, 3143, 3393, 3025,      4063, 6886, 7606, 8823, 8161.83, 8618.83, 7522.83, 6153.83, 4093.83,      3660.83, 3462.83, 3666.83, 4427.83, 6860.83, 8874.83, 9304.83,      9147.42, 9182.42, 7315.42, 6579.42, 4707.42, 4033.42, 3640.42,      3620.42, 4568.42, 7725.42, 8039.42, 9411.42, 10957.42, 8023.42,      9125.42, 6218.42, 5482.42, 4306.42, 3408.42, 3664.42, 4756.42,      5627.42, 9513.42, 12294.42, 11093.67, 10030.67, 10320.67, 6067.67,      5065.67, 5841.67, 4453.67, 3228.67, 5381.67, 7069.67, 8816.67,      12187.67, 12044.58, 9240.58, 8378.58, 6899.58, 6250.58, 4336.58,      2446.58, 4370.58, 3699.58, 7849.58, 9588.58, 10180.58, 10594.67,      9046.67, 9480.67, 7074.67, 4808.67, 3983.67, 4307.67, 4114.67,      5306.67, 7847.67, 10868.67, 11193.67, 10881.33, 10884.33, 9474.33,      6932.33, 4884.33, 4264.33, 4011.33, 4144.33, 4029.33, 7630.33,      10502.33, 11586.33, 11636.58, 10192.58, 9839.58, 6800.58, 4806.58,      4402.58, 4900.58, 3963.58, 4528.58, 7053.58, 8887.58, 10715.58,      11931.92, 10158.92, 9172.92, 7805.92, 5314.92, 3846.92, 5080.92,      4067.92, 5220.92, 6575.92, 10495.92, 12055.92, 12324.92, 9260.92,      8769.92, 8536.92, 4985.92, 3573.92, 3908.92, 4825.92, 5335.92,      7837.92, 9700.92, 12087.92, 12492.42, 11627.42, 9152.42, 6915.42,      5407.42, 4093.42, 3237.42, 3162.42, 4995.42, 7655.42, 9599.42,      11165.42, 12057.42, 10179.42, 10732.42, 6444.42, 5920.42, 4211.42,      3583.42, 3553.42, 4747.42, 7345.42, 9999.42, 11077.42, 11465.67,      11612.67, 10513.67, 6298.67, 5507.67, 4535.67, 4099.67, 4135.67,      4607.67, 6508.67, 9387.67, 12012.67, 15081.5, 12394.5, 11656.5,      6847.5, 5725.5, 3865.5, 3646.5, 4590.5, 4737.5, 6302.5, 8313.5,      9504.5, 10294.75, 9161.75, 8754.75, 5898.75, 5236.75, 4521.75,      4892.75, 4051.75, 6163.75, 7932.75, 10898.75, 11179.75, 10801.17,      9011.17, 10266.17, 8102.17, 5678.17, 5000.17, 4732.17, 4763.17,      6525.17, 7826.17, 9381.17, 10295.17, 10931.75, 9485.75, 8388.75,      6467.75, 4432.75, 4873.75, 5822.75, 4701.75, 4487.75, 7133.75,      8191.75, 10964.75, 12690.08, 10770.08, 9304.08, 7196.08, 6567.08,      5117.08, 4208.08, 4414.08, 5102.08, 7019.08, 9009.08, 12855.08,      11287.58, 10241.58, 9757.58, 6152.58, 6034.58, 4787.58, 4481.58,      4139.58, 4115.58, 6271.58, 8960.58, 9766.58, 10222.33, 12102.33,      7783.33, 7708.33, 5416.33, 4264.33, 4523.33, 3762.33, 5482.33,      6843.33, 8289.33, 10803.33, 11292.33, 10816.33, 13084.33, 8347.33,      5733.33, 4946.33, 4759.33, 5288.33, 6281.33, 5076.33, 8047.33,      10062.33), Weather = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 836.6924493,      864.3506976, 581.6976848, 517.7491028, 262.9073017, 125.8485446,      6.378750412, 47.44717309, 224.2027809, 564.7883676, 534.5681729,      804.9532423, 1008.232273, 661.7460732, 797.811819, 473.8225,      243.5474283, 196.7520236, 3.349023044, 66.17021767, 255.3949123,      273.7186901, 754.5074287, 1074.443123, 1138.053058, 1040.979912,      920.0242833, 500.2516119, 361.6047682, 148.0344631, 110.9923818,      59.62586443, 334.3492223, 460.5606125, 718.9706934, 1103.239148,      1135.812042, 680.8131205, 656.6975921, 603.0120248, 280.2865251,      113.6086427, 40.56542192, 13.99500372, 205.1454271, 537.9852507,      740.414947, 863.9388574, 840.1871901, 655.3627691, 689.800506,      450.142815, 214.0762313, 96.84368144, 94.40207546, 102.6299095,      219.5820493, 482.6409888, 849.0681028, 930.1270804, 836.8685145,      870.3884301, 666.1091718, 457.8937971, 223.2997102, 135.3552585,      29.08801136, 63.35303783, 56.79015527, 463.0870952, 753.1634369,      854.9329799, 926.6413953, 701.4458617, 685.8712135, 391.1133558,      172.7297975, 114.0125259, 118.8340228, 39.93146937, 195.3704427,      368.0819926, 608.7092419, 801.0635085, 938.8709503, 774.4973879,      747.6592081, 537.4781963, 192.260989, 188.0866469, 43.94658034,      36.07426405, 284.7563021, 271.1165662, 724.9144894, 998.2149659,      904.3262439, 618.0036623, 674.6240644, 513.5043282, 200.8462624,      67.53863402, 49.55548149, 13.91408483, 226.4356358, 512.7364786,      652.9808563, 977.2313968, 987.9207143, 980.2207144, 676.2825783,      481.8598638, 188.9155341, 18.81403734, 17.87685266, 22.92422788,      191.2999234, 628.7164688, 623.9494986, 877.1366222, 993.6142777,      780.0251185, 734.4778804, 433.5276702, 327.9133614, 143.3598735,      92.41551972, 43.18090657, 193.5962228, 418.4997331, 710.0634849,      918.6149901, 866.0983273, 943.8945311, 804.4710259, 451.2019972,      299.2649129, 130.4053496, 57.39215607, 109.9076214, 166.5745851,      363.4026804, 727.2489427, 942.2479917, 1140.126083, 920.1753988,      899.0587303, 506.487232, 258.806951, 121.5292418, 1.37176059,      108.8383996, 65.19810738, 305.5901759, 569.2363481, 731.8804303,      710.1873727, 683.3470574, 628.7776364, 348.321329, 233.8441238,      61.74325304, 79.11932747, 69.1111157, 257.976228, 518.1000209,      751.3500331, 884.8741716, 770.6473749, 705.6736025, 731.5307109,      524.8816439, 201.3512057, 79.29200986, 58.30966457, 46.91051312,      277.9127552, 466.5113866, 673.8802828, 908.1044645, 1101.706185,      858.9625893, 725.4699465, 315.9941158, 224.7973081, 163.5971151,      47.04241939, 25.25869003, 156.7798246, 519.0503168, 545.8445539,      963.5458059, 1193.839102, 917.7867611, 737.4186353, 480.8523752,      390.7623601, 116.7815813, 19.4804422, 76.28805273, 270.397472,      518.7399707, 697.1215068, 1193.349175, 933.5835268, 853.5872036,      713.9427207, 333.4000969, 215.9254222, 87.12224683, 96.06514093,      63.11633281, 149.1157803, 448.3832879, 704.5634245, 754.4911486,      872.9159352, 1050.3147, 598.0685504, 510.0915466, 214.3145038,      146.4499474, 60.86600378, 38.81855226, 206.9537231, 490.4757385,      697.0382613, 943.3598201, 985.8553275, 913.1275834, 995.459711,      503.2853052, 285.0997045, 134.6185906, 18.61433414, 41.02519091,      233.8134415, 386.077747, 696.2700349, 784.2664915), IndProduction = c(66.13584172,      66.75780888, 65.72119694, 65.2374447, 65.72119694, 59.50152531,      66.82691634, 66.13584172, 67.51799097, 67.24156112, 67.79442082,      68.2090656, 68.2090656, 68.27817306, 68.90014022, 68.13995813,      68.96924768, 69.79853724, 70.48961186, 68.69281783, 69.52210739,      70.35139694, 71.66443873, 68.96924768, 70.35139694, 70.62782679,      70.76604171, 71.11157903, 70.4205044, 70.76604171, 71.80265365,      70.55871932, 70.14407455, 70.35139694, 71.11157903, 70.21318201,      69.10746261, 70.48961186, 70.07496709, 71.04247156, 72.0790835,      70.76604171, 69.93675216, 71.45711634, 70.55871932, 71.66443873,      72.0790835, 72.14819096, 71.80265365, 72.35551335, 72.00997604,      72.49372828, 72.49372828, 73.80677007, 72.0790835, 75.05070439,      74.15230738, 74.42873723, 74.77427454, 76.43285365, 75.32713424,      75.05070439, 75.6035641, 77.88411036, 74.15230738, 77.12392827,      77.53857305, 77.33125066, 77.74589544, 77.95321782, 78.29875514,      79.61179693, 79.473582, 79.12804469, 79.88822678, 78.50607752,      80.51019394, 80.51019394, 81.27037603, 81.75412827, 82.51431036,      82.65252528, 82.65252528, 82.99806259, 85.35009976, 84.09758258,      83.5607895, 83.29239296, 82.57666885, 85.17116874, 83.29239296,      82.39773782, 81.41361718, 82.9345309, 83.91865155, 83.29239296,      84.45544463, 84.9027722, 83.47132398, 83.47132398, 82.75559988,      81.5925482, 80.60842756, 81.23468615, 80.16109999, 80.07163448,      78.9085828, 77.29820356, 76.2246174, 75.24049675, 75.68782432,      74.97210021, 75.06156573, 74.70370367, 73.36172097, 76.40354843,      75.77728983, 75.68782432, 75.59835881, 75.95622086, 74.8826347,      76.40354843, 76.31408291, 77.02980702, 77.02980702, 77.38766907,      77.20873805, 77.5666001, 78.10339318, 79.5348414, 79.80323794,      79.89270345, 77.38766907, 78.55072075, 78.55072075, 78.55072075,      79.26644486, 78.46125523, 78.19285869, 77.65606561, 78.99804831,      76.76141048, 77.65606561, 77.5666001, 77.5666001, 77.11927253,      77.74553113, 77.47713459, 78.37178972, 78.28232421, 77.47713459,      78.99804831, 79.44537588, 78.37178972, 79.35591037, 79.26644486,      78.46125523, 79.35591037, 79.89270345, 80.07163448, 79.17697934,      81.50308269, 82.12934128, 79.98216896, 80.69789307, 82.03987577,      82.30827231, 82.84506539, 83.47132398, 83.47132398, 84.2765136,      83.91865155, 84.09758258, 83.5607895, 85.26063425, 84.09758258,      83.20292744, 83.73972052, 82.66613436, 82.75559988, 83.91865155,      82.75559988, 83.11346193, 83.65025501, 84.00811706, 84.2765136,      84.63437566, 85.43956528, 85.35009976, 86.15528939, 86.33422041,      86.42368593, 85.11066398, 87.22334004, 87.32394366, 88.4305835,      90.54325956, 88.53118712, 90.04024145, 90.84507042, 90.94567404,      90.7444668, 91.34808853, 92.05231388, 91.75050302, 92.75653924,      91.24748491, 89.63782696, 90.64386318, 91.04627767, 88.32997988,      90.7444668, 89.53722334, 88.02816901, 86.92152918, 87.42454728,      87.92756539, 88.12877264, 88.32997988, 88.4305835, 87.52515091,      90.04024145, 88.32997988, 90.34205231, 89.23541247, 88.53118712,      90.14084507, 88.53118712, 89.23541247, 89.23541247, 89.03420523,      88.73239437, 88.02816901, 87.62575453, 89.73843058, 87.52515091,      87.22334004, 89.93963783, 90.7444668, 91.14688129, 90.44265594,      90.94567404, 90.34205231, 91.6498994, 92.65593561, 92.05231388,      93.15895372, 92.25352113, 91.95171026, 93.0583501, 91.54929577,      91.14688129, 93.9637827, 92.45472837, 93.56136821, 94.36619718,      93.15895372, 95.17102616, 96.57947686, 93.86317907, 96.47887324,      98.39034205, 97.18309859, 97.08249497, 97.58551308, 98.08853119,      97.18309859, 99.49698189, 100.9054326, 100.3018109, 102.1126761,      102.8169014, 102.7162978, 102.4144869, 104.2253521, 105.331992,      105.7344064, 106.2374245, 106.7404427, 105.8350101, 107.7464789,      107.6458753, 108.4507042, 108.6519115, 109.7585513, 109.8591549,      109.6579477, 110.8651911, 112.3742455, 112.0724346, 111.3682093,      111.6700201, 109.3561368, 110.5633803, 108.7525151, 110.9657948,      108.3501006, 105.9356137, 101.4084507, 97.48490946, 89.73843058,      86.92152918, 87.32394366, 84.7082495, 88.63179074, 89.3360161,      88.4305835, 89.63782696, 93.56136821, 91.54929577, 92.65593561,      92.65593561, 93.46076459, 92.85714286, 96.07645875, 97.48490946,      100.3018109, 100.2012072, 99.69818913, 101.1066398, 102.6156942,      104.527163, 104.527163, 105.8350101, 106.0362173, 107.444668,      108.249497, 108.5513078, 109.6579477, 108.0482897, 111.6700201,      110.4627767, 108.5513078, 109.5573441, 109.2555332, 107.5452716,      108.3501006, 108.4507042, 109.7585513, 107.444668, 109.3561368,      108.249497, 109.8591549, 109.1549296, 107.6458753, 105.9356137,      105.8350101, 106.9416499, 106.2374245, 107.0422535, 107.444668,      108.4507042, 107.0422535, 109.054326, 107.5452716, 109.8591549,      108.7525151, 107.5452716, 110.7645875, 110.1609658)), .Names = c("Demand",      "Weather", "IndProduction"), class = "data.frame", row.names = c(NA,      -360L))


title: How do you analyse an unbalanced dataset using dummy variables?

text: I have a returns (finance) dataset with 3 factors, each with 3 levels: industry (banking, retail, others), year (2010, 2011, 2012), and size (small, medium, large). Each factor has been coded as a group of dummies (e.g. IND_BANK, IND_RETAIL, IND_OTHERS). My sample is unbalanced, i.e. some cells do not have data at all. How would I go about coding it up as a linear regression with dummy variables and interactions? Here's a breakdown of a similar dataset. blob:http%3A//stats.stackexchange.com/b56e0691-09b7-41c8-a708-37e6c8c95944 I am thinking of modelling it as 3x3 dummies and 3x3x3=27 interactions; however, I clearly need to exclude some terms to prevent perfect collinearity. Also, I need to drop those interaction terms which have no cell data. Could anyone advise me on how to do this?


title: Estimating the dof for t test of Correlation Coefficients of smoothed time courses

text: I'm working with an array of size 265 x 5000. Each column is a time series and the rows are points in a series of volumes ( voxels in fMRI ). I preprocessed this data in several ways, most of which have some sort of filtering mechanism. My task is to compute the correlation coefficient between all voxels, compute a T test, find the p values and correct my coefficients for values above a threshold. Then there's more steps, but they are not relevant as far as I know. I know how to compute the T test of the corr coefs, but unclear on the DOF.               import numpy as np     # roi is my variable holding the 265 x 5000 data     corr_matrix = np.corrcoef(roi.T)# transpose of roi     corr.shape     >>>(5000,5000) # to be clear in case I messed up somewhere...      Thinking this might be helpful, I computed the rank of the corr array with:               np.linalg.matrix_rank(corr)     >>>225      The rank of the time series array:               np.linalg.matrix_rank(roi.T) # the transpose     >>>236      I don't really know what I'm doing here and hoping there's a clear way of estimating the dof of this array. Thank you! EDIT: I feel kinda dumb having computed rank. Neither matrix is a system of linear equations...so that probably made no sense...I think what I was after, is something akin to computing the eigen values of the cross correlation array. I've read in a genetic research paper dealing with large cross correlation matrices, that a fractional eigen value means there's dependence and integer's are independent. So that e.v. = 1.45 would mean 1 independent and 0.45 of another independent e.value. I haven't been able to confirm this. UPDATE: Maybe I wasn't wrong to compute the rank of the roi above. Since each col is a vector,if there are scalars a and b, which result in a*c_1 + b*c_2 = 0, then c_1 and c_2 are dependent. If roi has rank 236, then corr has dof = 236**2 ....?


title: When to use multiple models for prediction?

text: This is a fairly general question: I have typically found that using multiple different models outperforms one model when trying to predict a time series out of sample. Are there any good papers that demonstrate that the combination of models will outperform a single model? Are there any best-practices around combining multiple models? Some references:   * Hui Zoua, Yuhong Yang "Combining time series models for forecasting" International Journal of Forecasting 20 (2004) 69– 84


title: What is an instrumental variable?

text: Instrumental variables are becoming increasingly common in applied economics and statistics. For the uninitiated, can we have some non-technical answers to the following questions:   1. What is an instrumental variable?   2. When would one want to employ an instrumental variable?   3. How does one find or choose an instrumental variable?


title: What is difference-in-differences?

text: Difference in differences has long been popular as a non-experimental tool, especially in economics. Can somebody please provide a clear and non-technical answer to the following questions about difference-in-differences. What is a difference-in-difference estimator? Why is a difference-in- difference estimator any use? Can we actually trust difference-in-difference estimates?


title: Determine if an action increases the proportion of 1s in binary data with unknown population

text: We're trying to improve search results, and we're trying to determine if certain changes (adding a word to a synonym list, removing it from the query, etc) have a statistically significant improvement on the results for future searches. We currently have data for past queries and their matching "correct" results. We also have tests to sample from those past queries, perform them again, and look at the percentage of "good matches". A good match being one for which the known correct result is within the top recommendations by the search engine. For a concrete example: say we want to know if adding a synonym for show->display will improve results. We look through past queries and find all of those containing the word 'show', then we find the number of good matches. Then we add the synonym map, and test again for the same sample of queries. The results let us know how many queries improved, how many got worse, how many stayed the same, total number of samples, and total number of good matches (before and after the test). To get an idea of the numbers, we have about 260,000 past queries, and for most words there are less than 10,000 queries which contain those words (and thus we can test the entire sample). For this example the results are:               Sample size containing 'show': 5250      Good matches (base):           3212     Good matches (test):           3208     Improved:                      8     Got worse:                     12     Total change:                  -4      The question is which test(s) would be the most helpful in determining if this change is likely to represent the change in future unknown queries? (Using a 95% confidence level) Here's the things I've tried on R, and their results:   * Two sample T-test on both sets of data (represented by vectors of 1s and 0s) - not significant   * T-test on the difference of the sets - not significant   * Bootstrapping on the difference in mean of each set, and doing a T-test of that result - significant   * Bootstrapping each sample on their mean, then T-test two sample t-test of the results - significant   * Bootstrapping each sample on their mean, then t-test of the difference in the results - significant   * McNemar test - significant


title: What are some distributions over the probability simplex?

text: Let $\Delta_{K}$ be the probability simplex of dimension $K-1$, i.e. $x \in \Delta_{K}$ is such that $x_i \ge 0$ and $\sum_i x_i = 1$. What distributions which are frequently (or well-known, or defined in the past) over $\Delta_{K}$ exist? Clearly, there is the Dirichlet distribution and the logistic normal distribution. Are there any other distributions which come up naturally in this context?


title: How to test asymmetry?

text: I have a series and I need to check if the values are asymmetric or not (less or more) with respect to the mean. What test should I use?


title: What to read from the autocorrelation function of a time series?

text: Given a time series, one can estimate the autocorrelation-function and plot it, for example as seen below: ![The time series](http://i.stack.imgur.com/gp9Vp.png) ![ACF](http://i.stack.imgur.com/5sxja.png) What is it then possible to read about the time series, from this autocorrelation-function? Is it for example possible to reason about the stationarity of the time series? **Edited** : Here I have included the ACF of the differenced series with more lags ![ACF after differencing](http://i.stack.imgur.com/SKtnq.png)


title: Books about model selection in ecology

text: I am an ecology student and have to deal with 10 or 20 field variables, including species frequencies. I need to screen out what variables are most important in the occurrence of a bird species. What book would tell me the methods to do this?


title: How to simulate the success of an RF packet transfer knowing the packet size and BER

text: I'm trying to simulate a data protocol running over a 5 Mbps RF link. The datasheet for the RF chip gives a logarithmic BER vs SNR plot (BER := bit error rate; SNR := signal to noise ratio), but to simplify I am choosing to fix the BER to the same value that many of the RF characteristics are measured: BER = 10^-4. How do I determine the instantaneous success of a packet transfer of size N bits (N := { 36000, 18000, 9000, 4500})? I'm not looking for the probability of success (but that may be needed as an interim calculation). I need to be able to make a boolean decision about each packet as the simulation time ticks by. I also have a random number generator available as input and can keep a packet transfer history if necessary. I'm using the SimPy package for even simulation.


title: Is there a package for R that allows smoothing splines in GEE?

text: I run into a problem where I would like to build a GEE in R with cubic regression splines (or any other spline type) for a longitudinal data set and an urgent need for grouping and multiple autocorrelation structures. However, I did not find any package capable allowing for that. If anyone could give a suggestion for a suitable package, I would really appreciate that.


title: Half-normal distributed DV in generalized linear model

text: My dependent variable is, by origin, absolute residual* left after some regression; it is _distributed half-normally_. Now I plan to use Generalized linear model in SPSS (GENLIN) to regress it on some predictors (totally different from those which produced the residuals). What distribution type should I use for the DV? For continuous data, GENLIN offers Gamma, Inverse Gaussian and Tweedie distributions. Which to choose to model the half-normal? Or should I apply special transforms before? And what link function would be most appropriate? What can you advice on that? Thanks. *More precisely, I'm analysing positive residuals separately and negative residuals separately.


title: Machine learning with categorical data: Can I delete examples?

text: I have a machine learning project that uses a bunch of features to predict a class that has categorical values. The possible values are: 1, 2, 3, 4, and 5. I'm interested whether the class has the value 1 or not. Therefore I added a new binary variable that is 1 when the class is 1 and 0 otherwise. So far so good. The problem now is that my data set is quite small ~200 observations and the features for the values the class has are very close to each other. That is, it is very hard to distinguish between 1 and 2, but it would be way easier to distinguish for example between 1 and 5. Is there anything I can do do improve my results? For example, I'm really not interested whether when the instance I'm going to predict has the value 3. I just want a classifier that is good in recognizing 1 events.


title: Principal component regression and including variables not in the PCA

text: I performed a PCA on some data related to habitat (different types of vegetation cover, soil water etc.), but there were some variables that I did not include in the PCA because they were not strongly correlated with any of the other variables. I selected the first three principal components from the PCA whose eigenvalues were above 1. Now I want to run a multiple regression to analyse what variables will predict bird density. Is it valid to include my three principal components as well as the variables that I did not include in the PCA into this regression?


title: Strange values of Cohen's kappa

text: I have two raters who agree on 93% on the cases (two options: yes or no). However, when calculating cohens kappa through crosstabs in spss I get really strange outcomes like -0.42 with a sig. of 0.677. How can such a high agreement in percentage result in such a strange kappa? I don´t get it.


title: How to compare three different replacements for eggs in baking?

text: I am conducting a food science test to replace eggs with other alternatives in baking applications.   * Independent variable: three different formulas (substitution percentages)   * Dependent variables: texture profile, moisture, and sensory evaluation In terms of procedure, I would like to make a separation, my study will need to incorporate some texture analysis parameters. Let's say each muffin will need to be measured using a texture analyzer. Then I will need to set up a panel to rate overall acceptance, texture, and so on. Every subject will need to try the samples and asign a value (Likert scale). **What would be the best design in this case?**


title: 'Forward Stagewise' option in LARS algorithm

text: Can anyone help me understand the `forward stagewise` part in the LARS algorithm? I was reading the R code and could not figure out what is `updateR` and `downdateR` doing.


title: How to estimate unknown parameter via maximum likelihood?

text: For instance, we observed $n$ tossing of a biased coin with probability of heads being $\theta$. How can I calculate this parameter via maximum likelihood? How can I derive the log-likelihood formula and the correct maximum likelihoo estimate of $\theta$?


title: What is the relationship between $Y$ and $X$ in this plot?

text: What is the relationship between $Y$ and $X$ in the following plot? In my view there is negative linear relationship, But because we have a lot of outliers, the relationship is very weak. Am I right? I want to learn how can we explain scatterplots. ![enter image description here](http://i.stack.imgur.com/hVoLg.jpg)


title: Non-orthogonal technique analogous to PCA

text: Suppose I have a 2D point dataset and I want to detect the directions of all the local maxima's of variance in the data, for example: ![enter image description here](http://i.stack.imgur.com/OB1Mz.png) PCA does not help in this situation as it is an orthogonal decomposition and therefore cannot detect both the lines I indicated in blue, rather its output may look like the one shown by green lines. Please recommend any technique which might be suitable for this purpose. Thanks.


title: Why describe a sample as i.i.d.?

text: I’m hoping for some help understanding the concept of RVs with respect to their use in the theory of drawing inferences on a population from a sample. To draw inferences on a population using a sample it is said that the observations must be i.i.d. RV's. I'm considering the example of a weighted die. If you wanted to test if a die was weighted you could buy one and roll it 1,000 times. On each roll you could record the number generated in order to understand the probability distribution of the die. The outcome from rolling a die once can be represented as a RV that maps from the outcome of the roll to an integer in the set $\{1,2,3,4,5,6\}$. I know that a RV $X=f(x)$ is a real-valued function that maps from its domain to a subset of the real numbers so if I were to describe this experiment before it took place I may write down the vector $[X_{1},X_{2},…,X_{1000}]$ where $X_{i} \,\, \forall \,\, i\in[1,1000]$ are i.i.d. RVs. We can imagine rolling the die 1,000 times and then writing in the results for each experiment, yielding $[x_{1},x_{2},…,x_{1000}]$, the realizations of the RVs. Although $[X_{1},…,X_{1000}]$ are identical they seem to be regarded in texts as distinct RVs, and I'm wondering why? If we recognize that $X_{i}=f(x) \,\, \forall \,\, i$ (i.e. each $X_{i}$ is the exact same function), and the reason that $X_{i}$ need not equal $X_{j} $ for $i\neq j$ is because we fed a different input into $f(x)$, then isn’t it equally valid to say that $X_{1},X_{2},…,X_{1000}$ represent multiple observations on the exact same RV? Indeed, if you have two random variables $f(x)$ and $g(x)$, but $f(x)=g(x)$ and they have identical domains and identical ranges, then it seems confusing to argue that $f(x)$ and $g(x)$ are “different”? What seems to actually be happening is that you have one RV, $f(x)$, which is a way of describing the possible outcomes of rolling a dice, and on each roll you input a different element from the domain into the function and hence get a distinct output. So can anyone explain to me the intuition for describing this process as i.i.d. RV's, as opposed to different observations from the same RV?


title: data-transformation for multi-state analysis using traminer/mstate

text: I am trying to estimate a multi-state model and am following the setup in the book by Mills (2011). The research units can go through different states a, b, c, and d. The research units can move between any of the four states in either direction. My problem is how to set up my data to do so. I initially started out with the dataset in this format:               ID   year   state   v1   v2       1    2000   a       1    2001   a       1    2003   b       2    1990   a       2    1995   a        In Mills book on page 206, the suggested format is however:               ID   from   to   trans   Tstart   Tstop   time   status   v1   v2          where from to show all possible movements between states, trans a sequence of transition per id, Tstart and Tstop and Time are the time variables, status indicates whether one experienced a state or not (censored) and v1,v2 are covariates. I am doubting how I should set the data frame up, since in my case there are many movements between states in various directions. I started out with adding a column to my original data, so that I have data of type "SPELL" as it is referred to in the Traminer manual. However, by doing so and changing it to a say STS format, I also loose my covariates to which I need to hold on for the multi-state model. EDIT: In sum, how do I go from the format above to the format required by the mstate package?


title: Solution for error message in minitab

text: I have used minitab software to fit time series model and predict. For some models I got the following error message; "Fitted model may be nonstationary or noninvertible. Completion of computation impossible." Using another package the relevant models were fitted. What is the reason for that kind of error message and Is there solution to fix that error?


title: Forecasting volatility using GARCH

text: How do I forecast volatility using GARCH in STATA after estimating the conditional volatility?


title: Event classification in Higgs' Boson discovery

text: I am watching the seminar (around min 20) from CERN about Higgs' Boson discovery. They quickly go on event classification talking about _boosting decision trees_. I am not really aware on what they need to classify -I guess they assume the audience knows enough the project in depth. I wonder if someone knows the whole detection and classification process in subject. The problem statement is fundamental for me to understand it. _Update_ I am investigating on what features, classes and statistical measures have been used in the physics experiment and then feed to decision trees and multivariate analysis.


title: Exponentiate a quadratic form, which MVN?

text: I read that if you have a quadratic $f(x) = ax^2+bx+c$ and providing the leading coefficient $a < 0$, then $e^{f(x)}$ is the pdf of a normal distribution with mean $\mu = -\frac{b}{2a}$ and $\sigma^2 = -\frac{1}{2a}$. Does this generalise to $e^{f(\mathbf{x})}$ defining a multivariate normal distribution where $f(\mathbf{x}) = \mathbf{x}'\mathbf{A}\mathbf{x} + \mathbf{x}'\mathbf{b} + \mathbf{c}$. If so what are the conditions for this to be true and what is the mean and covariance? Thanks!


title: What is the $i$th sufficient statistic in the EM algorithm for Gaussian mixture models?

text: I am reading up on the EM algorithm for Gaussian Mixture Models, and there is consistent reference to the $i$th sufficient statistic. What is this, and why is it relevant to the algorithm?


title: Averaging several density functions

text: I have a data set grouped by a factor (like e.g. cyl in the built in mtcars in R). I plotted the estimated density funtions (desityplot) for each factor using the lattice package, e.g.               library("lattice")     densityplot( ~mtcars$mpg[mtcars$cyl==4] )     densityplot( ~mtcars$mpg[mtcars$cyl==6] )     densityplot( ~mtcars$mpg[mtcars$cyl==8] )      Now I would like to create a densityplot() for the whole dataset but the estimation should disregard the different number of datapoints available for each factor, e.g.               densityplot( ~mtcars$mpg )      is not what I want, because there are more data points with cyl=8 than data points with cyl=6. Is it possible to estimate the density separately for each factor and than calculate some sort of "average" density? If not, what is the appropriate method for my purpose?