title: PCA and FA example - calculation of communalities

text: I'm trying to understand how Principal Component Analysis and Factor Analysis work by implementing examples. Although I'm mainly using Python and Numpy here, this isn't Python-specific, as I'd like to know how to get the correct result generally speaking. There is an example in chapter 16 of Statistics in a Nutshell, but I can't understand how some of the resulting values are obtained. I'm not sure whether (a) there's something wrong with my implementation/calculations or (b) I'm misunderstanding the "jump" from PCA to FA. (The examples are on page 300-303. It wouldn't be appropriate to reproduce material from the book, but I hope it's reasonable to quote its short numerical example.) From the example table (psychometric test results) in the data matrix below:               import numpy as np          # (Table 16-1 in the book)     # Columns are:     # Reading, Music, Arith., Verbal, Sports, Spelling, Geometry     data = np.array(     [[  8,  9,  6,  8,  5,  9, 10],      [  5,  6,  5,  5,  6,  5,  5],      [  2,  3,  2,  6,  8,  6,  4],      [  8,  9, 10,  9,  8, 10,  6],      [ 10,  7,  1, 10,  5, 10,  2],      [  9,  8,  4,  9,  1,  7,  2],      [  3,  9, 10,  2,  6,  4,  9],      [  8, 10,  3,  8,  5,  7,  2],      [ 10,  9,  3, 10,  6, 10,  3],      [  7, 10,  1,  9,  6, 10,  2]])          # Transposing for np.corrcoef     data = data.T          corrmat = np.corrcoef(data)      The values calculated this way using NumPy for the correlation matrix match those in the book (Table 16-2 for those who have it). Here is the same input data for R:               data <- matrix(c(8, 9, 6, 8, 5, 9, 10, 5, 6, 5, 5, 6, 5, 5, 2, 3, 2, 6, 8, 6, 4, 8, 9, 10, 9, 8, 10, 6, 10, 7, 1, 10, 5, 10, 2, 9, 8, 4, 9, 1, 7, 2, 3, 9, 10, 2, 6, 4, 9, 8, 10, 3, 8, 5, 7, 2, 10, 9, 3, 10, 6, 10, 3, 7, 10, 1, 9, 6, 10, 2), ncol=7, byrow=T)     pc = prcomp(data, scale. = T)      NumPy's output for the correlation matrix is this:                 [[ 1.        ,  0.53484056, -0.25290194,  0.86021546, -0.46870501,          0.76225482, -0.38632342],        [ 0.53484056,  1.        ,  0.24875809,  0.26248971, -0.26308503,          0.38020761,  0.06879192],        [-0.25290194,  0.24875809,  1.        , -0.50102461,  0.20615027,         -0.30668865,  0.75803231],        [ 0.86021546,  0.26248971, -0.50102461,  1.        , -0.23649405,          0.89522479, -0.56880086],        [-0.46870501, -0.26308503,  0.20615027, -0.23649405,  1.        ,          0.05436758,  0.26604241],        [ 0.76225482,  0.38020761, -0.30668865,  0.89522479,  0.05436758,          1.        , -0.29078439],        [-0.38632342,  0.06879192,  0.75803231, -0.56880086,  0.26604241,         -0.29078439,  1.        ]]      I'm using the following to perform the PCA. The `data` matrix is turned into the `pca_output` matrix. The cummulative % also match the book's example (Table 16-4).               eigenvalues, eigenvectors = np.linalg.eig(corrmat)          # Order the eigenvalues by decreasing value     # (and then order eigenvectors).     evals_order = np.argsort(-eigenvalues)     eigenvalues = eigenvalues[evals_order]     eigenvectors = eigenvectors[:, evals_order]          # pca_output: columns are the principal     # components, after transposition (used here)     pca_output = np.dot(eigenvectors.T, data).T          cummulative_perc_variance = 100 * np.array(         [ eigenvalues[:i].sum()/eigenvalues.sum() for i in range(1, eigenvalues.shape[0] + 1)])      The first 3 eigenvectors I get are:                 [[-0.48310951, -0.2554222 , -0.08049741],        [-0.20665755, -0.60274754, -0.16463817],        [ 0.31165047, -0.56586657,  0.02618752],        [-0.5112334 , -0.00690029,  0.19728832],        [ 0.21569808,  0.04569528,  0.83434948],        [-0.43878073, -0.18261299,  0.4643702 ],        [ 0.35546879, -0.46450701,  0.12258562]]      I get the same results using R's `prcomp`:               > pc     Standard deviations:     [1] 1.8676341 1.2850439 1.0568965 0.6517665 0.4837756 0.2582021 0.1344177          Rotation:                 PC1          PC2         PC3        PC4        PC5         PC6         PC7     [1,] -0.4831095  0.255422196 -0.08049741  0.2324959 -0.1897195  0.76737470 -0.12638483     [2,] -0.2066575  0.602747536 -0.16463817 -0.6995409  0.2296628 -0.05498277  0.14750200     [3,]  0.3116505  0.565866565  0.02618752  0.1601325 -0.7081318 -0.22419271 -0.06802832     [4,] -0.5112334  0.006900293  0.19728832  0.2433956 -0.1357962 -0.28997811  0.73341720     [5,]  0.2156981 -0.045695275  0.83434948 -0.3293147 -0.1060965  0.33662557  0.14908327     [6,] -0.4387807  0.182612988  0.46437020  0.1413347  0.1812280 -0.38481239 -0.59798379     [7,]  0.3554688  0.464507012  0.12258562  0.4932350  0.5892964  0.11120226  0.19982732      The book then produces this table (Table 16-3: Communalities) and says: "The first step after computing PCA is to examine what proportion of variance is accounted for by the factor structure. This is done by examining the communalities [...]." $$ \begin{array}{r|r|r} & \mbox{Initial} & \mbox{Extraction} \\ \hline \mbox{Reading} & 1.0 & 0.929 \\ \mbox{Music} &1.0 & 0.779 \\ \mbox{Arith.} &1.0 & 0.868 \\ \mbox{Verbal} &1.0 & 0.955 \\ \mbox{Sports} &1.0 & 0.943 \\ \mbox{Spelling} &1.0 & 0.967 \\ \mbox{Geometry} &1.0 & 0.814 \\ \end{array} $$ It's not clear to me where these extracted communalities values come from. My understanding is that they are the sum of the squares of the loadings in a given row of the eigenvector matrix (over the first p columns where p is the number of selected components). According to the eigenvalues I get above, this would be $(-0.483)^2 + (-0.255)^2 + (-0.080)^2 = 0.305$ for "Reading" (first row) and 3 Components, and not $0.929$ (as in the table in the book). The book shows a table (Table 16-7) called "unrotated component matrix", as follows: $$ \begin{array}{r|r|r} & Comp 1 & Comp 2 & Comp 3 \\ \hline \mbox{Reading} & 0.902 & 0.328 & -0.085 \\ \mbox{Music} & 0.386 & ... & ... \\ \mbox{Arith.} & -0.582 & ... & ... \\ \mbox{Verbal} & 0.955 & ... & ... \\ \mbox{Sports} & -0.403 & ... & ... \\ \mbox{Spelling} & 0.819 & ... & ... \\ \mbox{Geometry} & -0.664 & ... & ... \\ \end{array} $$ This table does explain the communalities values ($(0.902)^2 + (0.328)^2 + (-0.085)^2 = 0.929$). My understanding was that the unrotated component matrix when doing FA was the same as the matrix of eigenvectors obtained for the PCA. From these examples, this doesn't seem to be the case. How should the "unrotated component matrix" be obtained and how does it differ from the PCA eigenvectors?


title: How to create a three dimensional matrix in matlab?

text: How to create three dimensional matrix? For example: A = matrix 4 * 5 (4 rows, 5 columns) B = matrix 4 * 10 (4 rows, 10 columns) How to turn these two matrices a three-dimensional matrix?


title: Simple relative frequency calculation, confusing wording

text: I am admittedly doing homework, but two problems are worded in such a way that is confusing me. Part (b) asks, "Determine the probability that an individual who has saved for retirement has saved less than $50,000. Use relative frequencies." Part (c) asks, "Determine the probability that a randomly chosen individual has saved less then $50,000 toward retirement." As far as I can tell, both parts are asking me to do the exact same thing. But that can't be right....


title: What are density and intensity of point pattern?

text: A simplified explanation with the focus on the following questions is being inquired. Appreciation goes in advance to whom provides **scientific** -- **simple** -- **practical** explanation.   Having an area **A** enveloping **N** number of points the density can be written as follows: $$D=N/A$$   **Question 1:**   What is the difference between **density** and **intensity**?   **Question 2:**   How to estimate density of a point pattern?   I'm aware of KDE methods but realized that there are serious difficulties in practice in regards to the bandwidth setting etc which are very effective in the estimation. To name I also tried **mass** function, 2D histogram, etc but am not satisfied yet with them. A robust, reliable and practical approach is being sought.   **Question 3:**   The following figure demonstrate another question. ![enter image description here](http://i.stack.imgur.com/T8N46.png)   Part **B** is the same as **A** removed some at the corner. Does it change the density of entire area? * * * to respected moderator: Is `question eligible for bounty in 2 days` a good mark? or! I'm receiving it immediately after my upload!


title: Normalizing weekly sales fluctuations in the overall market

text: I have a rudimentary question regarding how to "normalize" a set of time series data, and would appreciate your thoughts. To make it very simple, the hypothetical situation is as follows: Suppose we have the products A and B, both of which are of similar nature and functions and can be used as a substitute each other. We have weekly sales data of each product. There are 10 competing products out in the market, which are similar to our products A and B, and the overall weekly market data for those 10 competing products are also available. We want to know to what extent our products A and B tend to cannibalize each other, and would like to know the correlation between the two, based on their weekly sales data. However, the overall sales volume in the market (for our 2 products and the 10 competing products combined) tend to fluctuate from one week to another, and the fluctuations require us to adjust or normalize them before calculating the correlations between products A and B. Suppose that the weekly market fluctuations are indexed as follows:               Market:      Wk1  100 (=base)     Wk2  109      Wk3  103     Wk4  101     Wk5   96     Wk6  108     Wk7  117     Wk8  147     Wk9  135     Wk10 115       Then for the purpose of normalizing or adjusting weekly sales volume of the product A and B so that they would supposedly show sales numbers with overall market fluctuations being adjusted, we would calculate the adjusted weekly sales by this index. For instance, suppose the WK2 raw sales data for product A is 36, and then the adjusted WK2 sales for Product A is 33, which comes from 36 divided by 109 times 100. Does this make sense? or do we have better or more adequate ways to deal with the situations? Thank you for reading this for the rudimentary question and I would appreciate your thoughts. I would very much appreciate your thoughts.


title: I know the probability of an outcome in a sample, how to apply to population?

text: I admittedly am horrible at statistics. I know that the probability of an individual who has saved for retirement has saved less then 50,000 USD is 16.4% according to my sample data, but only 69% of workers say they have saved for retirement. So what is the probability for a random person selected from the entire population to have saved less than 50,000 USD? So that means 31% of the population admits to not saving anything (which is less then 50,000 USD), so the probability of a random person selected from the entire population having saved less than 50,000 USD for retirement is actually greater than 16.4% - but I don't know how to calculate what it is more precisely. This is homework, but I am so lost.... Any help is appreciated.


title: Simultaneous Non-parametric regression and Non-parametric density estimation

text: Given a collection of $(x,y)$ data, one might do a non/semiparametric regression of $y$ on $x$ to understand how to predict the $E(Y|X)$. Similarly, if one has enough data, it might be useful take a narrow vertical slice of the data (a small range of $x$), and try to understand the pdf of $Y$ at that point. However, I don't have very dense data, so I would like to combine these two techniques. What I would like to understand is, for any value of $x$, the full distribution of $Y|X$. I suppose that this really isn't even possible if the densities change rapidly as a function of $x$, but I'm willing to make some assumptions and take some short-cuts (I have a sense that the density is being 'stretched' as $x$ increases, and it should be unimodal though not symmetric) - I'm still trying to get a sense of what I need to assume. Any advice or references are appreciated. My only idea so far is to do it iteratively - regress, find the conditional means, estimate the densities, repeat. There seem to be some shortcomings, though.


title: Cross product term in linear regression why is it zero

text: ![regression](http://upload.wikimedia.org/math/8/d/e/8def8df3d62ea70cbf10a6e2b2b3127a.png) I'm trying to see why this is the case. I have taken the LHS and added $X\hat{\beta}$ and subtracted $X\hat{\beta}$. From this I can get both terms on the right but I end up with a cross product term $2(y-X\hat{\beta})^T(X\hat{\beta} - X\beta)$ that I can't prove equals to zero. I think there is some trick I'm missing, anybody remember it?


title: Chi-squared test of homogeneity

text: I am familiar with the chi-squared goodness of fit test and also the chi- squared test of independence. An example of the chi-squared goodness of fit is where we have frequencies for three different groups and want to know if the observed proportions differ from the expected proportions. An example of the chi-squared test of independence is where we have 2 or more categories measured within 2 or more other categories, and want to know which categories are associated. However I keep seeing a third type of chi-squared test called the 'test of homogeneity'. For example stattrek has a page on the chi-squared test of homogeneity: http://stattrek.com/chi-square-test/homogeneity.aspx Could somebody explain how the chi-squared test of homogeneity differs from the other two tests?


title: Introduction to support vector machines

text: I'm looking for a good introductory book/course on Support Vector Machines. My statistical background is almost non existent, so the more introductory the better.


title: Autocorrelation of a Bernoulli sequence

text: I have some doubts about the following exercise: Given the sequence of Bernoulli (0,1/2) independent random variables $x(n)$, $n\ge 1$, determine the average function and autocorrelation (time-delay) of the process: $z(n)=\sum_{m=1}^n y(m),\quad n\ge1$ with $y(m)=2x(m) -1 $ for significant values of time and delay. * * * I then calculated the average function: $E[z(n)]= n-1 $, Is that right? I have some problem for the autocorrelation: $E[z(n)z(n-k)]=E[\sum_{m=1}^n y(m) \sum_{m=1}^{n-k} y(m)]=...$ is it right to continue like this? i.e. What should I do with $z(n-k)$ in the summation?


title: How to calculate the percentage difference taking 95% CI into consideration?

text: I now have two averaged numbers, say $a$ and $b$. $a=1.1$ with 95% Confidence Interval of $\pm0.2$ and $b=1.0$ with 95% Confidence Interval of $\pm0.1$. Is it reasonable to calculate their percentage difference w.r.t $b$ as $$\%\ diff=\frac{1.1-1.0}{1.0}=10\%$$ From my understanding, due to the overlapping CI, we cannot realy say that $b$ and $a$ are different. But if I wish to quantify _how different they are_ , how may I do this reasonably with CI taken into consideration?


title: What plots should be used for diagnostics for linear mixed model?

text: Before fitting a linear mixed model, can any plots be used to show a random intercept/slope is justifiable in the model? I.e. these plots may indicate a different pattern for each individual over time. After fitting a linear mixed model, what diagnostics can be performed to see if any model modification is needed?


title: Error of data fitting by gamma.fit() in Python

text: I need to find gamma fit for data in Python 3.2.               param = gamma.fit(samp) // samp is a list of floating point numbers      I got error: `TypeError: cannot perform reduce with flexible type` Any help would be appreciated.


title: Can I use MLE here and what is the approximate SE?

text: I have got 10k data of the form `(y, x1, x2)` where `y` is the dependent variables and `x1, x2` are independent variables. I have an existing model `M`, where `M(x1, x2)` gives the pdf of `y`. So `M(x1,x2,y)` gives the density at `y` Now, I would like to improve the model by introducing some simple transformation of its argument, like: `M2(x1, x2) = M(a*x1 + b*x2 + e, c*x1 + d*x2 + f)` I picked the parameters `a ~ f` based on max loglik (`sum log M2(x1, x2, y)`) of the total 10k data set. Is this method of using MLE valid? (I don't need of cross validation at this point) it looks like nonlinear regression, seems a bit different. What's the distribution of `a^ ~ f^` approximately? Can I find a confidence interval of each parameter so that I know if they differ from 0 or 1 so I can simplify the model? (given `M` is continuously differentiable and I have got Jacobian and Hessian of it. Also I am able to get Jacobian and Hessian of the loglik from `M2` w.r.t. `a ~ f`)


title: Fitting ARMA model with MATLAB R2012b

text: I want to fit an ARMA model on a time series (quarterly log returns of a 10 year bond) using MATLAB R2012b. This is part of an exercise. I have problems with the code and the interpretation of a result. This is the code I'm using:               % last element is the most recent one, first element is the oldest one     timeSeries = ...         [0.0327; 0.0386; 0.0375; 0.0147; 0.0417; -0.0319; -0.0278; ...         -0.0693; 0.0328; -0.0721; -0.0727; 0.0645; -0.0301; -0.0301; -0.0261; ...         0.0313; -0.0912; -0.0235; 0.0436; -0.0714; 0.0128; 0.0506; 0.1451; ...         0.0387; -0.0190; 0.0822; -0.0009; -0.0094; 0.0408; 0.0324; 0.0661; ...         0.1519; 0.0713; 0.2015; 0.0353; -0.0800; 0.0079; 0.0119; -0.1376; 0.0345; ...         0.0568; -0.3770; 0.0516; 0.0563; 0.0849; 0.0484; 0.0148; 0.0991; 0.0796; ...         0.0872; 0.0754; -0.0473; -0.3189; 0.0478; 0.0457; 0.0183; -0.0422; ...         -0.0528; 0.0632; -0.0287; -0.1751; 0.0238; 0.0733; -0.0177; 0.0933; ...         0.1358; -0.0631; -0.0542; -0.0166; 0.0194; -0.0868; 0.0742; 0.0260; ...         0.0122; 0.0634; 0.0277; 0.0171; 0.0509; 0.1551; 0.0875; 0.0535; 0.0058; ...         0.1711; 0.1028; -0.2364; 0.0991; -0.0061; 0.0768; -0.0444; 0.2470; ...         0.0935; -0.0983; -0.0428; -0.0852; -0.0977; 0.0094; -0.3146; 0.1685; ...         0.0379; -0.1846; -0.3769; 0.0353; -0.1330; 0.2245; 0.0255; 0.1483; ...         -0.0051; 0.0255; -0.0407; 0.0779; 0.0305; 0.0353; 0.0931; 0.0520; 0.1156; ...         -0.0712; 0.0579; 0.0879];          p = 1;     q = 1;          % first way: garchfit     spec = garchset('R', p, 'M', q, 'VarianceModel', 'Constant');     [coeff, coeffErrors, llf, innovation, sigma, summary] = garchfit(spec, timeSeries);     garchdisp(coeff, coeffErrors);          % second way: arima class     [model, varCov, llf2, info2] = estimate(arima(p,0,q), timeSeries);      I get the following results using garchfit (coefficients equal to 1 or -1 have an absolute value very close to 1, like 0.999999998):                                              Standard          T            Parameter       Value          Error       Statistic       -----------   -----------   ------------   -----------                C    0.0063979      0.022847         0.2800            AR(1)    -1             0.025306       -39.5159            MA(1)    1              0.040576        24.6451                K    0.010616       0.0010517       10.0941      These are the results using the arima class:                                             Standard          t           Parameter       Value          Error       Statistic      -----------   -----------   ------------   -----------      Constant     0.00381023     0.0235065       0.162092         AR{1}       0.638413       2.19145        0.29132         MA{1}      -0.617482       2.26535      -0.272577      Variance      0.0107181    0.00103327         10.373      **Question 1**. I believe that the two ways to fit the model should lead to the same model: an ARMA model can be represented as a GARCH model with constant variance (this is what I was taught at least, I'm no expert of GARCH models). Is it correct? **Question 2**. Assuming that both ways should lead to the same model, why do I get very different values for coefficients AR{1} and MA{1}? Is it just related to the optimization procedure used to compute them? **Question 3**. Assumption: the garchfit results are correct; in fact, these results are significant, while the ones computed with the arima class are not. Therefore, my univariate model becomes $$ y_{t} = 0.0063979 - y_{t-1} + \epsilon_{t-1} + \epsilon_t $$ where $\epsilon$ are residuals. Therefore, it is $$ y_{t} + y_{t-1} = 0.0063979 + \epsilon_{t-1} + \epsilon_t $$ Can you help me understand what does it mean? I think it should mean that the time series is equal to its residuals, therefore the time series is a white noise according to the theory of ARMA models. Is is correct? Thank you very much.


title: What does "static" in the context of a SEM mean?

text: I just stumbled upon a wiki article: http://en.wikipedia.org/wiki/Endogeneity_(economics)#Omitted_Variable In the section "simultaneity" theres a description of a model which says that both variables are codetermined ,i.e. each affecting the other, and it's a _static_ SEM The model is $y_i = \beta_1 x_i + \gamma_1 z_i + u_i$ $z_i = \beta_2 x_i + \gamma_2 y_i + v_i$ **My question is:** What does static in this context mean? My guess is that though it's a dynamic model we actually can't observe the dependency between $y_i$ and $z_i$. One would need at least one additional "shifter" in order to observe these dependencies.


title: When can we speak of collinearity

text: In linear models we need to check if a relationship exists among the explanatory variables. If they correlate too much then there is collinearity (i.e., the variables partly explain each other). I am currently just looking at the pairwise correlation between each of the explanatory variables. **Question 1:** What classifies as too much correlation? For example, is a Pearson correlation of 0.5 too much? **Question 2:** Can we fully determine whether there is collinearity between two variables based on the correlation coefficient or does it depend on other factors? **Question 3:** Does a graphical check of the scatterplot of the two variables add anything to what the correlation coefficient indicates?


title: Significant differences and related variables

text: I am reading a paper in which a statistical difference has been found between two treatment groups for a variable A (p=.05), but not for variables B (p=.06) and C(p=.06), even though A=B-C. Is this possible, and if so, how? The paper states that "All data were statistically analyzed using the Student's independent t-test provided the data were normally distributed. Otherwise, the Wilcoxon-two-sample test was used."


title: Matrix of correlations between the residuals in a SEM fitted by R systemfit function

text: Assume we would fit this model with `systemfit()` from the systemfit R package. $y = \beta_1 x + \gamma_1 z + u$ $z = \beta_2 w + \gamma_2 y + v$ Now I can get an estimate of $Cor(u,v)$. I know that one cause of endogenity is that there is a significant correlation between $u$ and $v$. But does this mean that I should look at correlation matrix provided by the `systemfit()` output in order to conclude if there is evidence that this correlation might be the reason for the endogenity? Further if I control for the endogenity and the correlation between $u$ and $v$ decreases does this mean that some omitted variable is now controlled for such that the remaining error is just the residual noise?


title: Alternative Egger's test, without using standard error

text: A simplified regression equation $ES=\frac{a+b}{n_1+n_2}$ has been suggested as an alternative to Egger's regression equation $\frac{ES}{SE}=\frac{a+b}{SE}$, where ES=Effect Size, $n_1$=sample size of the patients, $n_2$=sample size of the controls, SE=Standard Error. This alternative test, that was presented by Peters et al. in their 2006 paper in JAMA, is supposed to be better than Egger's test when the ES is the lnOR. This alternative test could also be valuable in cases Standard Error (SE) cannot be calculated, as SE is not taking part in the equation. Could this alternative Egger's test be used with the other types of Effect Size? When the ES is the SMD? When the ES is the RR? When the ES is the Pearson's correlation coefficient?


title: Several questions about conditional probability

text: Suppose I have 2 variables $A$:   $P(A) =$ 0.01   $P( \lnot A) =$ 0.99 And $B$ that depends on $A$: $P(B|A) =$ 0.05   $P( \lnot B|A) =$ 0.95   $P(B| \lnot A) =$ 0.01   $P( \lnot B| \lnot A) =$ 0.99 Applying: $$P(B)=\sum_{A}^{ } P(B|A)P(A)$$ we get $P(B)=(0.01)(0.05)+(0.99)(0.01)=0.0104$ Ok, my question is the following: If I set the probability of $P(B)=1$   How do I get the values of $P(A)$? As $B$ depends on $A$, how are all the probabilities affected?   How to compute   $P(A)$ ?   $P(B|A)$?


title: How to compare the effectiveness of medical diagnostic techniques?

text: A pathologist friend came to me for help with the following question for a research project. The goal is to compare the effectiveness of three different diagnostic techniques. The data set is as follows: there are 50 different specimens, each specimen was evaluated by 4 pathologists, and 3 different instruments (ie 600 total diagnoses). Each case has a possible diagnosis of positive or negative, and the true results are known, as they have been independently determined. The success rate depends on both the quality of the instrument and the skill of the pathologist, and we cannot assume that the four pathologists have the same proficiency. Finally, even though each person measured the same specimen 3 times, they can be treated as independent measurements. What are the appropriate tests for comparing effectiveness among the instruments? Thanks. ADDED: Lots of good info in the answers, thanks to both. Any thoughts on how ROC and randomized block compare/contrast? I'm not sure I've digested it enough to know which method is "better". Since the results need to be communicated to a certain audience, it probably depends on which is more widely used among that audience.


title: Data preparation for regression

text: I am trying to predict real estate sales prices.   * In my dataset there are independent variables that are both nominal and numeric (square meters, prices etc.)    * Before feeding the data to any regression algorithm I'd like to preprocess it correctly (binning, normalizing mean / std deviation, discretization etc.)   * I am overwhelmed by the many methods listed in various textbooks and try to find out what works well in practice Although the most reasonable answer to this question is probably 'it depends', could you maybe give me some rules of thumb / war stories / general advice?   * How do you usually preprocess data for regression?    * What methods do you usually apply?   * Which regression algorithms need a special treatment? As my tools I am using weka and R. Many Thanks!


title: Which pseudo-$R^2$ measure is the one to report for logistic regression (Cox & Snell or Nagelkerke)?

text: I have a `SPSS` Output for a logistic regression. This output reports two measure for the model fit, `Cox & Snell` and `Nagelkerke`. So as a rule of thumb, which of these R² measures would you report as the model fit? Or, which of these fit indices is the one that is usually reported in journals? * * * Some Background: The regression tries to predict the presence or absence of some bird (capercaillie) from some environmental variables (e.g., steepness, vegetation cover, ...). Unfortunately, the bird did not appear very often (35 hits to 468 misses) so the regression performs rather poorly. Cox & Snell is .09, Nagelkerke, .23.   The subject is environmental sciences or ecology.


title: Metaregression with Age as a Covariate

text: Imagine that you have access to a collection of studies in which an age range, a sample size and an "effect" (say, cigarettes per day) is reported. What is the best way to construct a meta-regression of cigarettes per day by age? The hard part here is to correctly interpolate across these studies with different sample sizes and different age ranges. First, the distribution of ages won't be reported. Second, I only have the effect mean and variance for the entire study. In the limit, it seems reasonable that enough studies with enough age ranges will give me an appropriate estimate of the relationship between age and cigarettes per day. Do you know if this sort of analysis has been done before, and if the appropriate statistical tests have been derived?


title: The case for non-parametric models

text: I am relatively new to non-parametric models, such as the Hierarchical Dirichlet Process, Pitman-Yor processes, etc. but I have read that non- parametric methods are _data-driven_ , and that they can be more _robust_ than their parametric counter-parts. I am looking for evidence and an exposition of why this is the case. Perhaps more concisely, what are the benefits and drawbacks of non-parametric models (and what evidence there is of this). More informally, can it be argued that there is a reasonable trend to move towards non-parametric models in the machine learning / statistics literature?


title: Automated parameter selection for a GARCH model, in a similar manner to the forecast package

text: I was wondering: is there are a package in R for automated GARCH model selection? I'm thinking of something like what the forecast package does for ARIMA models. If I implement this myself, would it be appropriate to just do a grid search over the possible parameters for the GARCH and ARIMA parts of the model (using the rugarch package), and select the one with the lowest AIC (or BIC)?


title: explanatory variables may bias predictions

text: I' m asking this question out of sheer curiosity, my teacher was not able to explain it. If I'm using logistic regression with categorical variables they are coded like {1,2,3}. I guess it wouldn't change my results if I used {4,5,6}. But what if linearity of coding wasn't kept? (say {4,10,99})? What I deal with is just the way of coding factor variables but is it possible that the statistical inference could be "distorted" this way? Or, worst case scenario, would I be able to draw nonsense conclusions from categorical data just because they are coded in some way?


title: Initial conditions of experimental units in Completely Randomized Design

text: Since I had no answer to the second part of this question I would like to ask it again providing some more infos I gathered in the meantime. **Context**   We are designing a complete randomized trial in a forest. The plots (experimental units) exhibit natural variability in the response variables. Since we will set up 3 replicates for each of the 3 treatments we only have 9 plots and, depending on the random configuration of assignment of treatments to them, the initial conditions, as far as the response variables are concerned, can be significantly different.   In an effort to guarantee similar initial conditions, we have computed ANOVAs for all the possible initial configurations and selected a subset of only those whose $F$ statistic $\approx1$, thus equally partitioning total variance between and within groups. In our ED, for $F\approx1$, $p\approx0.42$: well above the conventional limit $\alpha=0.05$ needed for the rejection of null hyphotesis of no difference between treatments.   The actual configuration would then be randomly selected from the subset of $F\approx 1$ configurations. **Question**   The $F\approx1$ is an arbitrary threshold but one that makes much more sense than, for example, minimizing $F$, a strategy that would maximize variance within groups resembling a Randomized Complete Block design, without blocks!   Are there specific experiences in assuring similar initial conditions between experimental units? Do you reckon the $F\approx1$ _criterium_ is feasible? Thank you.


title: Modelling a mixed model in JAGS/BUGS

text: I am currently in the process of implementing a model for soccer result prediction in JAGS. Actually, I have implemented several, but I have reached my most difficult challenge yet: A model described by Rue & Salvesen in their paper "Prediction and retrospective analysis of soccer matches in a league". Their model uses a mixed model for truncating a Poisson distribution conditioned on attack / defense strength after 5 goals. They have also adapted a law from Dixon & Coles (1997) to increase the probability of 0-0 and 1-1 results in low-scoring games. My problem is as follows, I am trying to implement the mixed-model: $$ \pi_{g1}(x_{A,B},y_{A,B}|\lambda_{A,B}^{(x)},\lambda_{A,B}^{(y)}) = \kappa(x_{A,B},y_{A,B}|\lambda_{A,B}^{(x)},\lambda_{A,B}^{(y)})Po(x_{A,B}|\lambda_{A,B}^{(x)})Po(y_{A,B}|\lambda_{A,B}^{(y)}) $$ Where $x_{A,B}$ denotes the number of goals scored by the home team in the game between teams A and B, and $log(\lambda_{A,B}^{(x)})$ denotes the teams' strength. I have tried to implement these two laws in JAGS by using the zeros- ones trick, but with no luck so far (`error: illegal parent values`). My JAGS model so far:               data {         C <- 10000              for(i in 1:noGames) {             zeros[i] <- 0         }              homeGoalAvg <- 0.395         awayGoalAvg <- 0.098              rho <- 0.1     }          model {              ### Time model - Brownian motion         tau ~ dgamma(10, 0.1)         precision ~ dgamma(0.1, 1)              for(t in 1:noTeams) {             attack[t, 1] ~ dnorm(0, precision)             defence[t, 1] ~ dnorm(0, precision)                  for(s in 2:noTimeslices) {                 attack[t, s] ~ dnorm(attack[t, (s-1)], (tau * precision) /                                               (abs(days[t,s]-days[t,s-1])))                 defence[t, s] ~ dnorm(defence[t, (s-1)], (tau * precision) /                                                (abs(days[t,s]-days[t,s-1])))             }         }              ### Goal model         gamma ~ dunif(0, 0.1)              for(i in 1:noGames) {                  delta[i]            <-  (                                     attack[team[i, 1], timeslice[i, 1]] +                                      defence[team[i, 1], timeslice[i, 1]] -                                     attack[team[i, 2], timeslice[i, 2]] -                                      defence[team[i, 2], timeslice[i, 2]]                                 ) / 2                  log(homeLambda[i])  <-  (                                         homeGoalAvg +                                          (                                             attack[team[i, 1], timeslice[i, 1]] -                                              defence[team[i, 2], timeslice[i, 2]] -                                             gamma * delta[i]                                         )                                     )                  log(awayLambda[i])  <-  (                                         awayGoalAvg +                                          (                                             attack[team[i, 2], timeslice[i, 2]] -                                              defence[team[i, 1], timeslice[i, 1]] +                                             gamma * delta[i]                                         )                                     )                  goalsScored[i, 1] ~ dpois( homeLambda[i] )             goalsScored[i, 2] ~ dpois( awayLambda[i] )                  is0X[i] <- ifelse(goalsScored[i, 1]==0, 1, 0)             isX0[i] <- ifelse(goalsScored[i, 2]==0, 1, 0)             is1X[i] <- ifelse(goalsScored[i, 1]==1, 1, 0)             isX1[i] <- ifelse(goalsScored[i, 2]==1, 1, 0)             is00[i] <- is0X[i] * isX0[i]             is01[i] <- is0X[i] * isX1[i]             is10[i] <- is1X[i] * isX0[i]             is11[i] <- is1X[i] * isX1[i]                  kappa[i] <- (                             is00[i] * ( 1 + (homeLambda[i] * awayLambda[i] * rho) ) +                              is01[i] * ( 1 - (homeLambda[i] * rho                ) ) +                              is10[i] * ( 1 - (awayLambda[i] * rho                ) ) +                              is11[i] * ( 1 + rho                                     ) +                              1 -       ( is00[i] + is01[i] + is10[i] + is11[i]     )                         )                  # This does not work!             zeros[i] ~ dpois(-log(kappa[i]) + C)         }          }


title: Odds Ratio Between Overlapping Groups

text: Is it inappropriate to use odds ratio for groups that are subset of each other? For example, comparing GPA of female students to GPA of all the students. (This is just an example, not what I am really comparing) This is a part of a larger project, and my main interest in doing this is to show that the result will be different because of the groups used. And it's based on earlier literatures that I felt the need to show the distinction between the two groups before I proceed with my analysis. In my opinion, the odds ratio is just a descriptive statistics to compare odds between two groups and I do not really see why it would be wrong. However, Stata did not allow me to do it and so it leads to this discussion. Main question: why must the groups be dichotomous?


title: Anomaly prediction confidence for frequentist vs bayesian parameter inference

text: I am comparing the behavior of some implementations of Bayesian and frequentist approaches to parametric anomaly detection and currently trying to figure out the differences when the sample set is really small. The Bayesian case looks fairly simple - it is possible to estimate the probability of a given item being from the same distribution as the data obtained and cut it at certain threshold. However, the frequentist case is not so obvious to me. I have a feeling that the frequentist approach would be to refuse to answer whether an item is anomalous or not until we get enough data to reject a null- hypothesis that it is a normal item. For example, how would a frequentist statistician approach the following problem: > Given the data set [0, 0.5, -0.5] and knowing that the data comes from the > normal distribution, is it possible to estimate whether any of the following > elements [10, 5, 2] are anomalous with the significance level of 0.05 using > the frequentist approach? What happens step by step when the training data > set grows bigger, from 3 say to 100 items - at which point will it be > possible to evaluate the above given test items?


title: Data-driven removal of extreme outliers with Naive Bayes or similar technique

text: I have a large set of location data from a social network and would like to conduct a mobility study with it. For each _object_ , I have up to several thousand _locations_ where this object posted from. I quickly summarized properties of objects, which are   * location count (simply the number of locations I have of this object)   * distance traveled (all the distances between successive locations summed up)   * radius of gyration (standard deviation from the mean location, which is also termed the center of gravity)   * average velocity (self explaining) Now, in order to study mobility, I really depend of having only humanoid objects in the dataset. The problem is that, on the particular social network, there are also _non-humanoid objects_ contained (e.g. a radio station broadcasting always from the same location or an earthquake monitoring service broadcasting from the whole world). Thus, in the end, I have several classes of objects in my dataset, e.g. humans, stationary bots, dynamic bots, etc. My strong assumption is that these different classes show very different values in the feature space, i.e. an earthquake monitoring service should have high location count, high distance traveled, high radius of gyration, etc. Based on this assumption I wanted to use a simple supervised classifier (such as Naive Bayes) to classify my whole dataset so that I can **remove non- humanoid objects**. The training set would encompass about 200 objects which are gained by looking at the respective social network profiles of the objects. ![enter image description here](http://i.stack.imgur.com/gTy3Y.png) **N~2000** But as you can see in the above figure, the distribution of the above- mentioned variables are extremely long-tailed, i.e. there are some extreme outliers but they are sparse, and the distributions are certainly not normal, which would be required for normal NB. When I log10 transform the variables, some seem to approch a normal distribution (namely location count and distance traveled), but radius of gyration and average velocity are still not gaussian: ![enter image description here](http://i.stack.imgur.com/Tb4uo.png) **N~2000** Thus the apparent problem is that there hardly will be any outliers in a training set of just 10^2 objects (and acquiring a larger training set is not easily feasible). Thus the _priori_ probability of an outlier class will be very small.   * would this even make sense with this dataset? i.e. is there a classifier that can deal with these kinds of data or are the proportions of classes to unequal?    * One practical requirement would be that the technique is already implemented in some software (preferably R) because the removal of these outliers is just a small step in the huge process of my work...


title: Comparing impact of variable in a category for a given sample data

text: I am a software engineer and I have been asked to look at some performance issues that we are having with some batch-jobs on a system I am working on. The batch jobs generate data for different categories (A-Z), however the mix of categories varies depending on the time of day. The job is run 5 times a day. The total amount of data per job is approximately 2,000 entries, with a different mixture of categories (e.g. the 9am job has roughly 200 category B's, 1200 category A's and 600 category C's). I also have recorded the total runtime of the job, however I need a way of estimating/approximating the time it takes an individual data point for each category to run. I have data that shows the category (A-Z) for each run of the job (1-5) for the past 30 days. From using pivot tables I am pretty sure that a particular category affects the performance of the 9am job (the other jobs have no entries for this category and run well), however I need a way of "proving" this. I am unsure about how I should go about this. As I said, I can use a pivot chart to show that there are a lot of category B entries in the 9am job, however I would prefer some mathematically accurate way of supporting this hypothesis. Thankyou for your help. Aidos


title: Measure of variation around Median

text: Are there terms for measures of variation around the median? I know that we use "variance" for means of data but how can variance be defined for the median? Is there a similar statistic for median?


title: What is the RMSE normalized by the mean observed value called?

text: I have been using the `Root Mean Squared Error` (RMSE) to measure the accuracy of values predicted using a model. I understand that the value returned is using the units of my measures (rather than a percentage). However, I would like to quote my values as a percentage. The approach that I have taken is to normalize the `RMSE` by the mean value of my observations. Is there a term for `RMSE/mean` ?


title: How to do time series ( longitudinal) clustering based entirely on Shape of the curves?

text: I have a longitudinal (panel) dataset for investment growth for 120 countries covering the time from 1960-2008. Essentially it's viewed as 120 time series. What I am interested in is to group countries based on their shape of their growth curves over time. Thus whether they share similar Shape of their curves are the only criteria I need for grouping those countries. I have tried KmL package (K-means for Longitudinal Data), but it seems that (please correct me if I am wrong) this methodology produces the result that group countries exhibiting similar (investment growth) mean value (or magnitude), not exactly according to the similar shape. For example, KmL tends to group countries with high investment growth, median average investment growth, low investment growth, etc. The countries within those groups may have very different shape of curves over time. What I am looking for is regardless of the absolute value of investment growth. As long as the two countries exhibit similar pattern of their growth over time curve, they should be grouped together in one group. Could anyone tell me a way to implement this clustering? I have noticed from previous posts that cointegration test may work. Any suggestions will be greatly appreciated!


title: nested design after anova tukey

text: so I analyzed a nested design (split-plo, 2x4 factorial) with ANOVA and it indicated that there are differences in the levels of factor b. I went on with t-tests comparing the means depending on factor b. If the assumptions for t-tests are positive can I just compare 2 groups without considering the nested design? Or do I have to do ANOVAs with factor b being restricted to 2 groups each time? Did both and got different results. Thx B.


title: How does word co-occurrence let two words from the same topic link together in LDA topic model?

text: The essence of LDA is word co-occurrence. I want to know why? and does word co-occurrence mean two words appear together in a certain document or just appear in the documents collection together? Thank you all in advance.


title: How to estimate IBD (identity by descent) between two individuals mathematically

text: I work with data containing SNP genotypes for thousands of individuals. One of the task I had was to perform IBD between the individuals to find out the individual relatedness. I perform the IBD analysis with PLINK (http://pngu.mgh.harvard.edu/~purcell/plink/ibdibs.shtml). However, I always wondered how it works mathematically in a concrete way. I tried to find out on internet, but all I get is theoretically informations but not how you practically estimate the IBD parameters. So, maybe someone here can explain it to me? For example, let say I want to find out the relatedness between individual 1 and 2 with SNP genotypes:               Ind1: A/G A/A G/C T/G T/T     Ind2: A/G A/G G/G T/G A/T      How do I estimate the IBD given those genotypes between Ind1 and Ind2? Thanks!


title: Jeffreys prior for linear regression model

text: Consider the linear regression model $${\bf y} = {\bf X}\beta + {\bf e},$$ where ${\bf y}$ is an $n\times 1$ vector, $\beta$ is a $p\times 1$ vector, ${\bf e}$ is an $n\times 1$ vector. Assume also that $e_j\stackrel{ind.}{\sim} N(0,\sigma)$. What is the Jeffreys prior of the parameters $(\beta,\sigma)$? I am basically looking for a reference where I can find this.


title: Using deep learning for time series prediction

text: I'm new in area of deep learning and for me first step was to read interesting articles from deeplearning.net site. In papers about deep learning, Hinton and others mostly talk about applying it to image problems. Can someone try to answer me can it be applied to problem of predicting time series values (financial, internet traffic,...) and what are important things that I should focus if it is possible?


title: Understanding intercept in simple linear regression and why one variable is a predictor and the other is an outcome variable

text: I'm a new statistics student :) I have some questions about linear regression, i'm using R to do some tests. I have two simple lists, like:               > a <- c(1,2,3,4)     > b <- c(5,6,7,8)      then I do               > model <- lm(a ~ b)      the result (coefficient) is:               > coef(model)     (Intercept)           b               -4           1       the "relationship" is perfect (1). But i didn't understand the intercept value, why is it -4? Then if i change my R test with:               > b <- c(5,6,7,80)     > a <- c(1,2,3,40)     > model <- lm(a ~ b)     > coef(model)     (Intercept)           b       -1.0868825   0.5137503      > model <- lm(b ~ a)     > coef(model)     (Intercept)           a         2.125346    1.945622       obviously **lm(a ~ b) != lm(b ~ a)** But, what are there criteria for selecting which list to put in the first place?


title: Struggling to understand the relationship between two definitions of entropy

text: I am familiar with entropy as defined from information theory, which is simply $ H(X) = \sum_{i=1}^{n} p(x_{i}) \log p(x_{i}) $, for a discrete set of probabilities. In term of clustering, the entropy can be used to measure the quality of a clustering solution. It is defined as, $ E(S_{r}) = - \frac{1}{\log{q}} \sum_{i=1}^{q} \frac{n_{r}^{i}}{n_{r}} \log \frac{n_{r}^{i}}{n_{r}} $, where $S_{r}$ is a cluster, $q$ is the number of classes in the data set, $n_{r}^{i}$ is the number of documents of the $i$'th class that were assigned to the $r$'th cluster and $n_{r}$ is the size of the cluster. The definition is exactly the same in Empirical and Theoretical Comparisons of Selected Criterion Functions for Document Clustering. I understand the idea behind why it is used, although I don't understand how we got from the first definition of entropy to the second equation.


title: Does an r value of 0.60 mean that an increase in one variable is 60% likely to mean an increase in the other variable?

text: If one has an r value of 0.60, can one state that an increase in one variable is 60% likely to mean an increase in the other variable?


title: How to understand effect of RBF SVM

text: How can I understand what the RBF Kernel in SVM does? I mean I understand the maths, but is there a way to get a feeling when this kernel will be useful? Would results from kNN be related to SVM/RBF since the RBF contains vector distances? Is there a way to get a feeling for the polynomial kernel? I know the higher the dimension, the wigglier it is. But I'd like to get an intuition what the kernels do rather than try all possible kernels and pick the most successful.


title: Interpolation in multivariate time series

text: I have a problem in multivariate time series. The data consist of three time series related to foreign trade. Although my client is still doing research and attempting to find monthly data for all three series, she may only have quarterly data for some series. She has data for 10 years. I am considering interpolating monthly data from the quarterly data, and was wondering if this is a good approach, and what any drawbacks there are. Thanks


title: Error with Syntax for GLM with inverse.gaussian family in R

text: In my data set, I have 3 covariates ($V_1, V_2, V_3$) and response ($V_4$), generated by a program specifically written to generate random variables using an Inverse Gaussian distribution, in which $ \sigma^2 = 1 $ and $ \mu = \frac{1}{(\beta_0 + \beta_1 \cdot V_1 + \beta_2 \cdot V_2 + \beta_3 \cdot V_3)^\frac{1}{2}} $, where $ \boldsymbol \beta $ is the coefficient vector. When I use the command `fit <\- glm(V4 ~ V1 + V2 + V3, family = inverse.gaussian)`, R gives 2 errors:               Error: no valid set of coefficients has been found: please supply starting value     In addition: Warning message:     In sqrt(eta) : NaNs produced      The first one seems to be something specific to the Inverse Gaussian family for GLM, as I haven't gotten that error when I use the same command for other families. I don't know why the second is thrown, as all such $ \eta = \boldsymbol \beta^T {\bf V} $ lie in the interval $ (0, 8) $. If I manually use `start = c(1, 1, 1, 1)`, then it proceeds without errors but has the solution $ \hat{\boldsymbol \beta} = [2.491 \ 0.156 \ 0.572 \ 1.170 ]$, which is definitely not right as I used $ \boldsymbol \beta = [0.759 \ 0.825 \ 0.916 \ 0.972] $ to generate the data. The data set contains 1,000 tuples, so the randomness of the generation shouldn't allow the solution for $ \hat{\boldsymbol \beta} $ to vary by that much. When I have done this for other families, the difference bewteen any one element of $ \hat{\boldsymbol \beta} $ and $ \boldsymbol \beta $ hasn't been more than $ 0.05 $ for a data set of 1,000 tuples. What should I do to correct this and properly estimate $ \hat{\boldsymbol \beta }$? For those curious, I am running a GLM on a data set where I already know $ \boldsymbol \beta $ for the purpose of testing a program.


title: Generalization of cumulative probability models for ordinal Y

text: There are many models in existence for ordinal $Y$, for example the proportional odds ordinal logistic model, the continuation ratio model, and the cumulative probit model. The first and third of these are in a class of models specifying Prob$(Y \leq j | X) = F(\alpha_{j} + X\beta)$ where $F$ is a non-decreasing function (e.g., logistic or normal CDF). This can be generalized by having a family of functions (e.g., a 3-parameter family) and estimating the parameters of this family in addition to $\alpha, \beta$. Or perhaps there is a way to do this non-parametrically in $F$. Has anyone seen any methods developed along these lines? Among other things, this would be a competitor of quantile regression.


title: What is sequential change point detection?

text: Could someone point me to an introductory-level explanation (think Wikipedia entry) of what sequential change point detection is?


title: Is this a nested model and can I use the likelihood ratio test?

text: I have a data set, total_data and I applied a model to it. For instance, the model has one parameter $\beta$, and I calculated the log-likelihood of the fitted model (using maximum likelihood method). In the mean time, I have a stratification variable that could divide total_data into two subsets: high_risk_data and low_risk_data. Now, I could apply the same model to each subset, obtaining $\beta_1$ and $\beta_2$. I could sum up the log-likelihood from each fitting and obtain some kind of overall log-likelihood for the total_data under model 2. Model 1: parameter $\beta$ fitted on total_data;   Model 2: parameters $\beta_1$ and $\beta_2$ fitted on each subset of data. Can I consider them as nested models, since I have one fewer parameter in model 1? Can I apply the likelihood ratio test to select a better model? I know a more standard way to obtain something similar to model2: introduce a second stratification variable-group and apply the model ($\beta$,$\beta_{group}$) to the total_data. The likelihood ratio test would be to test against H0:$\beta_{group}=0$. Is this the same as the method I described above?


title: KMO value of 0.06

text: I have a dataset that contains data from 307 subjects and nine variables for each subject. I would like run a PCA. My problem is that I get a KMO value of 0.06. My question is whether this is due to the fact that within every subject the variables add up to 1. Am I right with this assumption? and is there another way to test whether its ok to run PCA on my data? Thank you in advance Matti


title: How do I compare multiple arbitrary predictions for a given data set?

text: Background: I am developing a Python Statistics Framework, not because the ones out there are bad but because it will help me learn Python and Statistics. I have taken AP Stats, and read scattered books and articles beyond. I am more than willing to read up on whatever shiny technique will solve my problem: The issue is that I don't know the name of said technique yet. Problem: Given two or more things, each of which takes in an X value and returns the probability of that result(Within a certain fixed range, so .00 to .01, .01 to .02 would each be separate blocks catching all the x values within that range and returning .005 and .015 respectively), and a data set, quantitatively figure out which function best matches the data. Doing the reverse(taking in a probability and returning an X value) would be a bonus. Idea: Be able to compare Logistic Regression, a Data Tree, and "If yes within the past 3 years then .8 else .01" style predictions. Is there a sane way to do this? Thank you.


title: Compare measure vs. true value: Methods and Visualization

text: I have a simulation routine that produces measurements of a continuous variable and I also know the true value. What is the optimal method to graphically depict, analyse and present the error or measurement? Any references, links, ideas, hints to the right direction are more than welcomed.


title: Why doesn't the exponential smoothing forecast package in R provide confidence intervals for the fitted values?

text: The upper and lower prediction intervals for the forecast periods are provided by the forecast() function. However, neither prediction or confidence intervals seem to be available for the fitted values within the range of the actual data. Why is this?


title: Using regression splines for values outside of the calibration range

text: in my question on a load forecast model using temperature data as covariates I was advised to use regression splines. This really seems to be a/the solution. Now I face the following problem: if I calibrate my model on winter data (for technical reasons calibration can not be done on a daily basis, rather every second month) and slowly spring arises I will have temperature data outside of the calibration set. Are there good techniques to make the regression spline fit robust for values outside of the calibration range? Any experiences or references? Thanks!


title: Observed info matrix via Hessian

text: In some resources, I saw that the observed information matrix is the negative of the expected value of the Hessian matrix. However, in some other resources I saw that it is just the negative of the Hessian. Also, in R this is calculated by the inverse of the Hessian matrix (no negative). I am extremely confused. Can you please explain this?


title: Variables lack correlation, but have pattern

text: Below is the graph of two variables, X and Y, each representing count data. N=348. Note the scales of the axes:   ![http://i.imgur.com/tNGyTX5.jpg](http://i.stack.imgur.com/ih1yU.jpg) Y is very approximately lognormal, but X has no decent fit (including Poisson, negative binomial, lognormal and gamma of the log transform).   Spearman coefficient between X and Y is close to 0, and p-value to reject no correlation is very high. From the plot, there appears to be no combinations of extreme values of both x and y. When I log transform both X and Y, the following plot results:   ![enter image description here](http://i.stack.imgur.com/DVVHU.jpg)   Clearly any appearance of pattern has disappeared. My questions are:   * Why is there a lack of combinations of "extreme" values on the linear scale, but not on the log scale?   * Is there any significance to the lack of combination of extreme values on the linear scale, and is there anyway to investigate further?  The purpose of this study is exploratory.


title: Some Basic things we need to do when we are doing text classification

text: I am working on a project where I have to do multi-label text classification. I want to understand that whether my approach is correct or I am missing something. I am using R to do it. >   1. Clean the text >   2. Create a corpus. While creating the corpus I am removing the sparse > term with sparse value set to .86 >   3. Create a DTM from this corpus and attach the label to it. >   4. Divide the DTM into training and test set >   5. Build the model using the train set >   6. Test the model using the test set. >   7. Do statistical analysis to prove the model > So my questions are >   1. Do we have to remove sparse term while creating the corpus and if so > does the sparse value is too high and is there any range of sparse value to > be set in text classification or it depends on the result we get. >   2. Is it a thumb rule that the number of rows should be more than the > square of the number of features. And also does number of features depend on > the number of unique classes you have >   3. Can we do feature selection once we do sparse term removal. >   4. Which one gives good result dimension reduction or feature selection. >


title: How to interpret the results from polr() function in R, ordinal logistic regression

text: I would like some help on how to interpret my results for the `polr()` function in `R`- for my ordinal logistic regression. It is particularly hard to interpret results as there are no p values. I have tried the the standard model-fitting functions, including predict, summary and do not understand the results it gives. Can some explain stepwise way to interpret results after completing a `polr()` model? I have read multiple times the post but feel that it does not clearly explain how to interpret results: How to understand output from R's polr function (ordered logistic regression)?


title: Fitting Crossover design as nested ANOVA model - Residuals issue

text: I am re-doing an analysis done in a paper _Experimental Designs in measuring the effectiveness of promotions_ where they have fitted a nested ANOVA model to a crossover designed sample ![enter image description here](http://i.stack.imgur.com/O372d.png) In the paper they have fitted period and treatment (ABC) nested in replicates. `R` model looks like               aov(sales ~ replicate + store + replicate/treatment + replicate/period)      But when I plotted the residuals in each treatment the box plot showed exact same variance across levels (too good to be true.) ![enter image description here](http://i.stack.imgur.com/0AW1W.jpg) I tried to observe the residuals, and notices residual sum in replicates = 0 within replicates sum of residuals for trt =0, period =0 and store is 0. What am I missing?


title: Separation in logistic regression in a complex survey?

text: Firth's penalized maximum likelihood estimates, exact logistic regression and Bayesian logistic regression (e.g. bayesglm) can account for separation in logistic regression. But how to account for probability weights and complex survey designs in these methods (e.g. in R or Stata)?


title: Interpretation of intercept of a regression line in time series data

text: Does the intercept value of a regression equation have meaning in a time series dataset? Suppose I have a dataset: the intercept is 27.512, but we are 95 percent sure that the intercept is between -34 and 89.074. Does the intercept value bear meaning? **Statistics using SPSS:** ![enter image description here](http://i.stack.imgur.com/afEFA.jpg) ![enter image description here](http://i.stack.imgur.com/IY3QN.jpg)


title: What is the proper test for inter-rater reliability (three raters, ordinal data)

text: I need help for choosing a proper test for inter-rater reliability. I have three raters (experts) who have analysed the quality of images (10 images from different location). The data is ordinal (image quality is assessed numerically from 1 to 5).


title: Can SVM be used for dimensionality reduction from n dimensions to n-1?

text: Consider a case when I have data that has $n$ dimensions and I want to transform this data into $n-1$ dimensions. For this problem in PCA we project all data points to a $n-1$ dimension plane preserving the variance. If I draw a hyperplane using SVM (that will be a $n-1$ dimension plane) and project all data points to that plane, how will it differ from the PCA plane?


title: Why do two equivalent approaches to selecting numbers that sum to 100 produce different results

text: This question is inspired by two seemingly equivalent attempts to solve a question asked yesterday: In R, how to sample from the output of combn(a,b) if the "a choose b" is too large? Essentially the problem is to choose $10$ integers that sum to $100$, where each number must be $\ge 0$. The approach suggested in the original question is to sample $9$ numbers in the range $1$ to $109$, _without replacement_ , and treat these as cut points cutting the range $0$ to $110$. The differences between cutpoints (including $0$ and $110$) give 10 numbers each of which is $\ge 1$. Finally, we subtract $1$ from each number to get $10$ numbers that are $\ge 0$ and add up to $100$. An equivalent -- and simpler -- approach, I thought, would be to sample $9$ numbers in the range $0$ to $100$, _with replacement_. The differences between cutpoints (including $0$ and $100$) directly give 10 numbers that are $\ge 0$ and add up to $100$. Implementing these approaches in R shows, however, that the approaches are **not** equivalent.   * `y` is the first approach, sampling $9$ numbers from $1\ldots109$ without replacement and adjusting.   * `x` is the second approach, sampling $9$ numbers from $0\ldots100$ with replacement. From the histograms below, it is apparent that $0$ is sampled a lot more often in the first compared to the second approach. ![enter image description here](http://i.stack.imgur.com/PymJJ.png) Assuming there isn't simply a bug in my second approach, my questions are:   * Why the difference when it seems like intuitively the two approaches should be equivalent?   * Which approach is correct?  R code below for reference               set.seed(99)          N=1e5     y <- sapply( 1:N, function(i) {         diff( c(0,                  sort(sample(1:109, 9)),                 110) ) - 1     } )          x <- sapply( 1:N, function(i) {         diff( c(0,                 sort(sample(0:100, 9, replace=TRUE)),                 100) )     } )          p1 <- hist(y, breaks=-0.5:100.5)     p2 <- hist(x, breaks=-0.5:100.5)     par(mfrow=c(3,1))     plot(p1, col=rgb(0,0,1))     plot(p2, col=rgb(1,1,0, 0.75))     # http://stackoverflow.com/questions/3541713/how-to-plot-two-histograms-together-in-r     plot(p1, col=rgb(0,0,1), main="x superimposed on y")     plot(p2, col=rgb(1,1,0, 0.75), add=T)


title: How to diagonalize a covariance matrix into Identity matrix

text: Given a symmetric, positive definite, covariance matrix $S$, how can I calculate a matrix $C$ such that $CSC^T$ is the Identity matrix $I$?


title: How to calculate the specific Standard Error relevant for a specific point estimate within a linear regression?

text: My objective is to actually get to the Standard Error for every point estimate of a linear regression. I originally thought that the Standard Error of a linear regression model was constant. And, that it was simply equal to the square root of the Mean Squared Deviation of the residual. That's why the model's Standard Error is also referred to as the RMSE. I now understand that when you construct a Confidence Interval around the trendline of estimates... at each point estimate, the applicable Standard Error is different from the overall model's Standard Error or RMSE as described above. The first part of my question, is how do you calculate this specific Standard Error at a specific point estimate? Could one provide a link to a numerical example to facilitate my interpretation of the formula? I understand that such Standard Errors have a predictable pattern. They are much lower than the model's overall Standard Error for values fairly close to the Mean value. And, they are larger from the same model's Standard Error as you go further away from the Mean value. As a result, Confidence Intervals are not parallel to the regression trendline. Instead, they curve like an hourglass shape. My second part to this overall question is: How come the resulting hourglass shape of the resulting Confidence Interval as depicted does not break the linear regression assumption that the variance of residuals remain constant across observations (the heteroskedasticity thing)?


title: Skewed data with many zeros: Simple and multiple regression

text: I read some research papers about Lachenbruch methods- 2 parts model. I am a clinician cum researcher and don't have much knowledge/experience in statistics. We used a scale (St. George's Respiratory Questionnaire) in a study (in general village population in India), 0 indicates no problem and 100 indicates severity. The dataset has a large number of zeros (indicates no problem), and is positively skewed. We have to do simple and multiple regression analyses to determine the associated factors. We use Stata, but we are not experts. Could you please let us know the process of analysing this data? Could you please help us with the Stata commands?


title: What classifier to choose?

text: Hi all I am doing a supervised classification of a binary labelled data where I have a mixture of categorical and continuous attributes. I have in total 6 attributes and nearly 6500 instances. I was wondering if there is an expected performance of a classifier based on (bias/variance) that might influence my decision to pick one classifier over another. I have missing data, noisy labels and outliers in my data. The list of classifiers I am trying to choose from are Logistic Regression, kNN, decision trees and ensemble learners. I tried weka and so far k-nn has given the best AUC after 10-fold cross- validation and testing on a new dataset. I am trying to find given the strengths and pitfalls of these classifiers if one is better over the other and why? Highly appreciate any help.


title: How do I calculate random baseline?

text: I am a bit confused as to how to calculate random baseline. If I understand correctly the random baseline is calculated by adding up the squared probabilities of all the classes. The random baseline classifier thus picks a class at random, instead of choosing the most frequent one. I have 7 classes, each with # of items and a total of X. How do I find the probabilities?


title: Similarity between [x,y,z] and [0,0,0]

text: I have two vectors $v_1$, $v_2$. Each can have arbitrary values and I measure their correlation. $v_1$ and $v_2$ can also be different in length and I zero- pad the shorter to macht the longer vector length. Shift invariance of correlation allows to find the the shorter within the longer vector. Sometimes either one can have all attributes set to 0. If $v_1 = [1,2,3]$ and $v_2 = [0,0,0]$, correlation returns NAN. Mathematically it makes sense, but not in my context. $[0,0,0]$ is a valid entry and I still want to detect the similarity of $v_1$ and $v_2$. How can I capture the correlation or similarity when $v_1$ or $v_2$ is 0?


title: Calculating chi-squared values for percentages?

text: I'm using scipy and I'd like to calculate the chi-squared value of a contingency table of percentages. This is my table, it's of relapse rates. I'd like to know if there are values that are unexpected, i.e. groups where relapse rates are particularly high:                       18-25    25-34    35-44    ...     Men     37%      36%      64%      ...     Women   24%      25%      32%      ...      The underlying data looks like this:                       18-25           25-34           35-44          ...     Men     667 of 1802     759 of 2108     1073 of 1677   ...      Should I just use the raw values, so have a contingency table like this, and run a chi-squared test on the raw values?                       18-25        25-34        35-44      ...     Men     667          759          1073       ...      That doesn't seem quite right, because it doesn't capture the relative underlying size of each group. I have been Googling, but haven't been able to find an explanation I understand of what I should do. How should I find unexpected values in data like this?


title: Gini index vs entropy

text: If I have a discrete probability distribution $p$ with $K$ classes   * Gini index = $\sum_{K}$$p_k$(1-$p_k$)   * Entropy = -$\sum_{K}$$p_k$log$p_k$ Per 'The Elements of Statistical Learning', `"Both the Gini index and cross- entropy are lower for the second split. For this reason, either the Gini index or cross-entropy should be used when growing the tree.` " I assumed two classes($K=2$) and probability values of $p_1$=0.1 and $p_2$=0.9. I got $Gini-index = 0.0081$ and $cross-entropy = 0.146$   * What conclusions can I make of this?   * Is there a formal proof that Gini index is a better approximation of misclassification error? What approach does one take for proving it?


title: Confidence interval vs Significance level question

text: Until now I thought that the significance level and the confidence interval were the complement of each other, i.e. when some question asked me to carry out a hypothesis test with a $p %$ significance level I would construct a $(1-p)$ confidence interval and see if it englobes or not the null hypothesis. However I was making a problem along the following lines that made me doubt some of my concepts: The problem gave me some random variable probability density function and some sample values and asked to carry out hypothesis test for $H_0:\theta=1$ and $H_a:\theta\neq1$ at the $5\%$ significance level. Since the alternative hypothesis is a composite one I used the likelihood ratio test together with the maximum likelihood to find a relation in the form $$\frac{L[H_0]}{max L[H_a]}<k$$ Since the pdf was a very unusual one I used Wilks' theorem and approximated it asymptotically by a $\chi^ 2$ distribution with 1 degree of freedom, so I could have a value of $k$ that represented the significance level I wanted. Until now everything went smoothly. However in the next part of the problem it asks me to outline briefly how I might obtain a $95\%$ confidence interval for $\theta$. Correct me if I am wrong but this wasn't what I made in the last example? Feel free to not only explain this specific problem in particular but to also give me a broader definition of the confidence interval and significance level. Sorry if I made any silly mistake. I started to study statistical inference in the last couple of weeks and I am still struggling even with its most simply definitions.


title: Whitening a dataset with fewer observations than variables

text: I have a k x n dataset where k equals the number of variables and n equals the number of observations per variable. I know these data are correlated and I would like to whiten them with the ordinary whitening transformation. Inconveniently, k outnumbers n by far, so that when estimating the k x k covariance matrix it will not be invertible. To get around this problem, I estimated the covariance matrix using optimal shrinkage, but the obtained covariance matrix is not suited for whitening anymore. I'd be grateful for ideas on how to whiten these data.


title: Visualising the variance

text: I work in a company that gathers work environment surveys. When a survey is done we create reports that are handed out to the managers in the company to show where they need to focus their efforts and such. In these reports we show a comparison of the variance of the calculated group compared to the variance of a reference population. We "normalize" this score on a 0-100 scale. (The in-code calculation gives another scale from -inf to +inf, but for simplicity we say that 0 is 50, and just chop off everything over 100 and under 0.) Our biggest problem is that our customers gets very conscious of the actual number. Even though we try to tell them that a high/low number isn't necessarily bad, it just shows how your groups variance compares to the reference populations variance. We are thinking of moving away from showing the number in the end-user report, and going for a visual representation of "high", "normal", "low" variance instead. But I can't figure out any good visual representation of it. It needs to be kind of neutral looking, but still show something... (.. i know ..) Does anyone have any suggestions on how this could be achieved? (Disclaimer: I am not a statistician, I'm a developer ;). I had one class of Statistics in my higher education, and that is over 5 years ago. So both the terms I use and my explanations could be totally meaningless.)


title: Does there exist a likelihood based experimental methodology book on undergraduate level for psychology / social science students?

text: I've been looking for an introductory book in statistics and experiment design for a methods course for psychology students. I think I've been looking through at least 30 books and all go through the same old motions. Fist hypothesis testing, then t-test, Anova, two-Way Anova and ordinary least square regression. And all of them insist on explaining these models by tedious calculations. A t-test requires you to calculate a t-statistic, an Anova requires you to sum squares of this and that, regression instead is done my minimizing the squared error. From my perspective I find it _much_ easier and intuitive to understand these types of models from a likelihood perspective. For example, a t-test is the model $y \sim \text{Normal}(\beta_0 + \beta_1x,\sigma^2)$ and an Anova is $y \sim \text{Normal}(\beta_0 + \beta_1x_1 + \beta_2x_2 + \ldots,\sigma^2)$ where $x, x_1, x_2$ are dummy coded variables. To fit such a model one would then use, for example, maximum likelihood which I find intuitive and easy to explain. This kind formulation also makes the connection between a t-test, an ANOVA and regression trivial to understand. So what I've been looking for an introductory book in statistics and experiment design that teaches statistics from this likelihood perspective. I would be very grateful for any suggestion!


title: Simple Mixed Model with 1 Fixed and 1 Random Effect

text: I have various datasets I need to analyse regarding soil properties, all in the same fashion, with one fixed effect (which is a position along a transect, indicating different land uses). Now my main level of replication is across different transects, which will obviously have some form of random variance associated with them, and so I want to account for this in my statistical analyses. So, in `lme4` I specified a mixed model to this specification               model <- lmer(variable.of.interest ~ transect.position + (1|transect))      Now, when I analyse the above model against a model without the transect position term, I get exactly the results I was expecting, and then plugging the above model into `anova()`, I get the F values, d.f. etc. that I need. However, I can't figure out how to say, for my overall report, that the random effect of transect does not make any difference to the overall analyses (i.e. I can't get a p-value, F value, d.f. etc.). Help?


title: Unexpected factor structure in different samples

text: I have administered a questionnaire to respondents asking about customers' adoption of technology. I want to do an exploratory principal components factor analysis in Stata to form constructs from a couple of questions (in line with the existing literature). The variables that I want to construct are perceived usefulness, perceived ease of use, perceived risk with product (PRP), etc. and I am comparing customers from two countries: The Netherlands (n=63) and Bulgaria (n=101). I designed items based on the research literature. So I expected that there would be no problems with the factor analysis; i.e., all items would have high factor loadings on a single factor. I have obtained cronbach's alpha which is 0.7 for all factors and 0.6 for usefulness for the Dutch customers, which is OK. The factor analysis is fine for most of the factors and most of items load most on one single factor as suggested by literature. The problem is that for PRP I got 4/5 factors for Bulgarians and 6/7 factors for Dutch which are actually different constructs compared to PRP for Bulgarians (13 factor loadings for PRP). I also have similar problem with Subjective norm (SN). My questions:   * When items do not fit into one factor is it possible to just take the means of such items? *   * Or is there something else I should do?    * Can you also provide me with some reference?


title: Should you report variance of asymmetrical data, such as ratios?

text: I have measured the time taken to solve a problem by algorithm $X$ and by algorithm $Y$. It takes a quite long time, so I have only 10 data for each algorithm: $$ X : ( x_1, x_2, \dots , x_{10}) \\ Y : ( y_1, y_2, \dots , y_{10}) $$ In my paper I reported the ratio $$ r = \frac{\sum_{k=1}^{10} x_k }{\sum_{k=1}^{10} y_k } $$ and provided a Wilcoxon signed rank test results. However the reviewer of my paper asks for **variance across the 10 ratios**. The distribution of a **ratio** is intuitively highly **asymmetrical around 1.** (you have only the interval $(0; 1)$ to capture the fact that algorithm $X$ is faster, but the entire $(1 ; \infty)$ to capture the fact that $Y$ is faster). So even a well estimated standard deviation can be of little use. This question is closely related to my previous one. The accepted answer proposes to **symmetrize** the data using the **logarithm** function and to construct **confidence** intervals for the ratio. However I have very many ratios (since I tested on many problems) and I am afraid that writing confidence intervals for each ratio would decrease the readability of the paper.   1. Should I decline to provide the variance of ratios?   2. Is there some other way to report the variance other than standard error?   3. If I decline, is there some literature which mentions that it is a bad idea to report the variance of ratios?


title: Does it work better to subtract mean of data in logistic regression?

text: I am using logistic regression to predict $X \rightarrow Y \in \{0, 1\}$ based on the featurization $\phi(X)$. The training objective function is \begin{equation} \mathcal{L} = \sum\limits_{i=1}^{n}{y_i \log \sigma(\theta^\top \phi(x_i)) + (1-y_i) \log (1-\sigma(\theta^\top \phi(x_i))}, \end{equation} where $\sigma$ is the sigmoid function. Since the objective is strongly convex, optimal solution can be computed via gradient ascent. However, I observe the following:   * The gradient is $\partial \mathcal{L} / \partial \theta = \sum_{i=1}^{n} (1-\sigma) \phi(x_i) y_i - \sigma \phi(x_i) (1-y_i)$.   * If the label set is unbalanced, e.g. 90% of $y_i = 1$, then most time the gradient is $\sum(1-\sigma) \phi(x_i)$.    * If the features $\phi$ are not zero mean, then most time $\theta$ gets positive updates.  However, if some feature $\phi_k(x)$ is negatively correlated with label, then the parameter $\theta_k$ should be negative. How could this happen?


title: Meanings of "1" on arrows on Confirmatory Factor Analysis diagram

text: I am just a novice in factor analysis. I am writing some presentation about factor analysis and I am mainly reading about exploratory factor analysis. Both EFA and CFA have diagrams. For EFA usually factor loadings are seen on the arrows. However, for CFA usually there are some "1" on some arrows, but not all. As I am a novice in this, I do not really understand yet what those mean. Can you give me some clues? For example, see here page 174, there is a picture : http://psychology.concordia.ca/fac/kline/Library/k13b.pdf


title: EM for Mixtures of Bernoulli (M-step)

text: When applying the M-step for a mixture of Bernoulli distributions, one of the parameters in our maximization is the Bernoulli parameter $\mu_{k}$, where $k$ is the index of the "mixture component", and $$ p(x|\mu_k) = \prod_{i=1}^D \mu_{ki}^{x_{i}}(1-\mu_{ki})^{(1-x_{i})}. $$ In our maximization with respect to this parameter, we get the following expression $$ \begin{align} \frac{\partial}{\partial \mu_{ki}}\mathbb{E}_{Z}[\ln p(X, Z | \mu, \pi) &= \sum_{n=1}^N \langle z_{nk} \rangle \left( \frac{x_{ni}}{\mu_{ki}} - \frac{1 - x_{ni}}{1 - \mu_{ki}} \right) \\ &= \frac{\sum_n \langle z_{nk} \rangle x_{ni} - \sum_n \langle z_{nk} \rangle \mu_{ki} }{\mu_{ki}(1-\mu_{ki})} \end{align} $$ where $$ \langle z_{nk} \rangle = p(z_{nk} | x_n, \mu_k, \pi_k) $$ Obviously, setting this to zero and solving for $\mu_{ki}$, we get the standard solution $$ \mu_{ki} = \frac{\sum_n \langle z_{nk} \rangle x_{ni}}{\sum_n \langle z_{nk} \rangle} $$ With that in mind, my question is as follows. Isn't there a constraint on $\mu$ such that $\sum_i \mu_{ki} = 1$? If so, then why is this not included in the maximization; i.e. why don't we formulate the Lagrangian which includes the term $ \lambda \left(\sum_i \mu_{ki} - 1\right) $?


title: method/regression to detect signal envelope

text: Are there methods to get the trend of a signal envelope? An example of signal here : ![enter image description here](http://i.stack.imgur.com/jolMw.png) The black line is the signal, the red lines are what my eye see as the envelope of that signal. I'd like to know what method could be use to get these envelopes and to quantify their relevance.


title: How to set up the initial data in the nls model

text: I am struggling with nls model given below. I have a blur idea how to set up the initial values in it. Shall I use "try and error" method or is there any other way to make it more systematic.               > model3<-nls(H2~H1*(A2^b1*(A1^b1*(b3+H1+sqrt((b3+H1)^2+2*H1*exp(b2)/A1^b1))+exp(b2))/     + A1^b1*(A2^b1*(b3+H1+sqrt((b3+H1)^2+2*H1*exp(b2)/A1^b1))+exp(b2))), start=c(b1=0.1, b2=2, b3=1), trace=T)      I get response               Error in nls(H2 ~ H1 * (A2^b1 * (A1^b1 * (b3 + H1 + sqrt((b3 + H1)^2 +  :        singular gradient      or               Error in nls(H2 ~ H1 * (A2^b1 * (A1^b1 * (b3 + H1 + sqrt((b3 + H1)^2 +  :        step factor 0.000488281 reduced below 'minFactor' of 0.000976562      depending on what values of parameters I use. Models formula look like this: ![enter image description here](http://i.stack.imgur.com/38MUg.jpg) Best regards


title: Removal of outliers

text: I work in a chemical laboratory and endeavour to estimate Measurement Uncertainty (MU) from sample retests. The procedure is well established. But what I an dealing with is a set of data where some duplicates are quite legitimate and reflect the underlying MU while others are not as the first test of the pair was, for some reason, incorrect. A double distribution, as it where. I tend to favour the method of using the '1.5 x IQR' method to eliminate outliers. I have found that often the trimmed data, when reassessed by the IQR method, has further outliers. The question is, in my circumstances, is it appropriate to apply this in a second round on the already trimmed data? Am I approaching the true underlying MU of the test or just fooling myself? (As I type this I have decided to create a dummy data set of a combined normal distribution and a skewed distribution and see what happens.)


title: Analytical formulation of a Hierarhical Bayes problem

text: In the free online book Bayesian Methods for Hackers, the last figure shows the estimation of the expected value of $\lambda$ for any given day: ![ente for image description here](http://i.stack.imgur.com/jkoAp.png) It looks like the author is calculating the expected value of the posterior for $\lambda$ for a given day. To do so, the author re-uses the MCMC samples from the full estimation of $\lambda_1$ and $\lambda_2$. But how does the **posterior of $\lambda$ for a given day** relate to $\lambda_1$ and $\lambda_2$? And why can one **re-use** the samples from the full MCMC calculation? In case it matters, below is the Python code in question that calculates the red line in the graph above.               figsize(12.5, 5)     # tau_samples, lambda_1_samples, lambda_2_samples contain     # N samples from the corresponding posterior distribution     N = tau_samples.shape[0]     expected_texts_per_day = np.zeros(n_count_data)     for day in range(0, n_count_data):         # ix is a bool index of all tau samples corresponding to         # the switchpoint occurring prior to value of 'day'         ix = day < tau_samples         # Each posterior sample corresponds to a value for tau.         # for each day, that value of tau indicates whether we're "before"         # (in the lambda1 "regime") or         #  "after" (in the lambda2 "regime") the switchpoint.         # by taking the posterior sample of lambda1/2 accordingly, we can average         # over all samples to get an expected value for lambda on that day.         # As explained, the "message count" random variable is Poisson distributed,         # and therefore lambda (the poisson parameter) is the expected value of         # "message count".         expected_texts_per_day[day] = (lambda_1_samples[ix].sum()                                        + lambda_2_samples[~ix].sum()) / N


title: How to perform genetic-algorithm variable selection in R for SVM input variables?

text: I'm using the kernlab package in R to build an SVM for classifying some data. The SVM is working nicely in that it provides 'predictions' of a decent accuracy, however my list of input variables is larger than I would like and I am unsure as to the relative importance of the different variables. I'd like to implement a genetic-algorithm to select the sub-set of input variables that produces the best-trained/fittest SVM. I'd like some help with choosing which R package to use when attempting this GA implementation (and possibly a brief psuedo-example). I've looked a most of the R GA/P packages out there (RGP, genalg, subselect, GALGO), but I'm struggling conceptually to see how I would pass in my ksvm function as part of the fitness function and input my variable array as the population pool...? Any help, thoughts, or nudges in the right direction gratefully received. Thanks code that solves this added below in a later EDIT               # Prediction function to be used for backtesting     pred1pd = function(t) {     print(t)     ##add section to select the best variable set from those available using GA       # evaluation function - selects the best indicators based on miminsied training error     mi.evaluate <- function(string=c()) {         tmp <- data[(t-lookback):t,-1]         x <- string         tmp <- tmp[,x==1]         tmp <- cbind(data[(t-lookback):t,1],tmp)         colnames(tmp)[1] <- "targets"         trainedmodel = ksvm(targets ~ ., data = tmp, type = ktype, kernel="rbfdot", kpar=list(sigma=0.1), C = C, prob.model = FALSE, cross = crossvalid)         result <- error(trainedmodel)         print(result)         }          ## monitor tge GA process     monitor <- function(obj) {     minEval = min(obj$evaluations);     plot(obj, type="hist");     }          ## pass out the GA results; size is set to be the number of potential indicators     gaResults <- rbga.bin(size=39, mutationChance=0.10, zeroToOneRatio=10, evalFunc=mi.evaluate, verbose=TRUE, monitorFunc=monitor, popSize=50, iters=3, elitism=10)          ## now need to pull out the best chromosome and rebuild the data frame based on these results so that we can train the model          bestChro <- gaResults$population[1,]     newData <- data[,-1]     newData <- newData[,bestChro==1]     newData <- cbind(data[,1],newData)     colnames(newData)[1] <- "targets"     print(colnames(newData))          # Train model using new data set     model = trainSVM(newData[(t-lookback):t, ], ktype, C, crossvalid)     # Prediction     pred = as.numeric(as.vector(predict(model, newData[t+1, -1], type="response")))     # Print for user inspection     print(pred)     }


title: Does my data come from a gamma or beta distribution?

text: I have data and I want to assertain whether it is beta or gamma distribution. Once I know what the distribution is, how do I find out what the parameters are? btw I am a complete novice


title: Parameter Search for SVM on the whole data

text: I am trying to implement SVM and i did my parameter selection(grid search) on the whole data and used the best values of C and gamma from that search to test on the testing data. Sometimes, the cross-validation accuracy, i get is 100%. Is this method wrong? What is the correct way of doing? How should i go about to do my parameter selection? Shld i divide the data(100%) into 10% for grid search and remaining 90% into training and testing? Need some guidance on it.. Thanks.


title: How to retrieve model-specific information when using the caret package?

text: CARET is a convenient tool that wraps dozens of modelling methods. CARET allows the user to access performance statistics such as RMSE, confusion matrices, etc. in a generic way without having to recall the usage of various modeling methods. However, some models have unique reporting outcomes. For example, with a boosted tree model one can compare the in-sample _vs._ out-sample error-rate as a function of the number of trees to test for overfitting. Naturally this report would not make sense with other methods that do not use trees. Is there a way I can use CARET to access these model-specific reports (maybe by accessing the modeling object generated directly)?


title: Residual analysis of cross-sectional time-series forecasts

text: I have forecasts and actuals for panel data (i.e. time-series cross-sectional data). The forecasts are already generated and provided by some source outside of R. I'd like to evaluate the quality of the forecasts. Are there standard tools in R that perform various diagnostics on the residuals? By diagnostics I mean tests such as:   * auto-correlation of residuals across the cross-section   * auto-correlation of residuals along the time series for a given member   * tests for fixed effects vs. random effects   * heteroskedasticity, etc. Or is the best way to perform these diagnostics to perhaps build a panel model using the forecast as the predictor in the panel model?


title: What to do when both dependent and some independent variables are repeated measures?

text: Good afternoon I am working on a study where the dependent variable (% of people in one ethnic group voting, essentially) and one of the independent variables (% of such people in the population) are measured annually for 10 years. There are other independent variables that are not time dependent. I could run a separate model for each year of the dependent variable, including only the earlier years' data for the independent variable, but was wondering if there was some more parsimonious method. I have access to R and SAS (If I've overlooked something obvious, my apologies, I have a bad cold and am not at my best) Peter


title: what should be the parametric form of the l2 regularization in a Bayesian setting?

text: In a Bayesian setting for parameter estimation, what should be the parametric form of the prior distribution in order to perform l2 regularization?


title: Can I have too much data?

text: I am training a large set of neural networks for a quite simple task. 10 of the networks have the same configuration, but have different amount of data. The 10 networks each have one hidden layer, with 2 neurons in it. The first network gets 1.000 training examples, the next 2.000 and the last one gets 10.000 training examples. The 1000 training examples for for the first network, is a subset of the 2000 training examples for the second network, which again is a subset of the 3000 training examples for the third network and so forth. I train my networks using the build in NN-Toolbox in MATLAB, where I use the Levenberg-Marquardt algorithm. When I train the networks, they all end up with a mean-squared error around 0.007 (which seems fair for my particular problem). The only one which differs, is the first one, which achieves a MSE of around 0.002. After I have trained my networks on the data, I test them on some test data which I did not use for training. Now the problem is, that the network with 1000 training examples seriously outperforms the other networks. From the MSE- score, this seems fair, but I do not understand why the networks with more available data, and even the data which network one uses, fails to learn the parameters better. Is there such a thing as too much data?


title: Does good calibration surpass roughly met assumptions and mediocre discrimination?

text: My question arises from my current task to develop a clinical prediction model using ordinal logistic regression (with rms), but applies to any kind of regression analysis.   1. The proportional odds assumption (difficult to assess...) is only roughly met. (I am aware of the continuation ratio approach).   2. Results of measures of discrimination (another difficult topic) are not overwhelming.   3. But the calibration plot is quite nicely. So, do I have to worry about 1 and 2 or does the very good calibration surpass basic assumptions and discrimination?


title: F ratio and model comparison method to compute confidence bounds

text: I am working with Fisher Ratio. I compute this ration from a function that was already implemented in the software I am working on. But I would like to be sure that it does what I think it does. Prism proposes a method for confidence intervals computation using model comparison. This uses F ratio(page 109). With dfd(=df2-df1)=1 and dfn(=df2)=99, I find a F of 1.08 approx with my code. Does it sound possible? I am using a library function computing the upper limit of integration of the f-density for the Cumulative Distribution Function F distribution. Inputs are confidence level and degrees of freedom for numerator and denominator. Is this going to give me F ratio? thanks and regards


title: Book recommendations for reviewing probability and statistics for engineering major

text: I'm an Industrial Engineering Major and I'm looking out for book(s) on Probability and Statistics (Combined would be preferred but separate would do as well) which :   * I, as a **self learner** , can work with.   * Is concise but covers : _Statistical methods for data analysis. Descriptive plots and statistics; observational studies and experiments; confidence interval estimation; hypothesis testing; regression and correlation._ in depth. **I have already studied them once, I just need to review them with greater depth.**   * Cater to engineers and not statisticians. I would prefer theory and application rather than theory and rigor.   * Would help me in my goal of studying Stochastic Models/Simulations/Financial Engineering (again, on my own).   * Has an international edition.   * Recommended/Not Required : Gives examples in R. I've been through the similar question but 2 answers are ones which recommend a book which the author himself isn't happy with and the accepted answer is a website. I've looked at Sheldon Ross, DeGroot, Rohatgi and many more but I'm still not sure if any of them give me what I want.