title: Restriction matrix for a VAR

text: In _New Introduction to Multiple Time Series Analysis_ by Luetkepohl (2005), section 5.2.1, it says that one can specify linear restraints for a VAR, $Y = \beta X + U$, in the form $$ \operatorname{vec}{\beta} = R\gamma + r $$ where $\beta$ is a $K(Kp+1)$-vector of parameters, $R$ is a $K(Kp+1) \times M$ matrix, $\gamma$ is an $M$-vector of "unrestricted parameters," and $r$ is a $K(Kp+1)$-vector of specified constants. He then derives an equivalence with the "standard" form of the restriction, $C\operatorname{vec}{\beta}=c$. Suppose I have a VAR like $$ \left(\matrix{w_t\\x_t\\y_t\\z_t}\right) = \left(\matrix{\beta^{ww}_{t-1}&\beta^{xw}_{t-1}&\beta^{yw}_{t-1}&\beta^{zw}_{t-1}\\\beta^{wx}_{t-1}&\beta^{xx}_{t-1}&\beta^{yx}_{t-1}&\beta^{zx}_{t-1}\\\beta^{wy}_{t-1}&\beta^{xy}_{t-1}&\beta^{yy}_{t-1}&\beta^{zy}_{t-1}\\\beta^{wz}_{t-1}&\beta^{xz}_{t-1}&\beta^{yz}_{t-1}&\beta^{zz}_{t-1}\\}\right) \left(\matrix{w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}}\right) + \left(\matrix{u^w_t\\u^x_t\\u^y_t\\u^z_t}\right). $$ How would I specify $R$, $\gamma$, and $r$, and/or the equivalent $C$ and $c$ (both is better), to obtain $$ \left(\matrix{w_t\\x_t\\y_t\\z_t}\right) = \left(\matrix{\beta^{ww}_{t-1}&\beta^{xw}_{t-1}&0&0\\\beta^{wx}_{t-1}&\beta^{xx}_{t-1}&0&0\\0&0&\beta^{yy}_{t-1}&\beta^{zy}_{t-1}\\0&0&\beta^{yz}_{t-1}&\beta^{zz}_{t-1}\\}\right) \left(\matrix{w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}\\w_{t-1}&x_{t-1}&y_{t-1}&z_{t-1}}\right) + \left(\matrix{u^w_t\\u^x_t\\u^y_t\\u^z_t}\right)? $$ And how would I extend it to the case with more than one lag, e.g. $$ \left[\matrix{\tilde{\beta}_{t-1} & \tilde{\beta}_{t-2}}\right] = \left(\matrix{\beta^{ww}_{t-1}&\beta^{xw}_{t-1}&0&0&\beta^{ww}_{t-2}&\beta^{xw}_{t-2}&0&0\\\beta^{wx}_{t-1}&\beta^{xx}_{t-1}&0&0&\beta^{xw}_{t-2}&\beta^{xx}_{t-2}&0&0\\0&0&\beta^{yy}_{t-1}&\beta^{zy}_{t-1}&0&0&\beta^{yy}_{t-2}&\beta^{zy}_{t-2}\\0&0&\beta^{yz}_{t-1}&\beta^{zz}_{t-1}&0&0&\beta^{zy}_{t-2}&\beta^{zz}_{t-2}\\}\right)? $$ In particular, I'm interested in testing $$ H_0: \beta = \beta_0 = \left[\matrix{\tilde{\beta}_{t-1} & 0}\right], \quad H_A: \beta = \beta_A = \left[\matrix{\tilde{\beta}_{t-1} & \tilde{\beta}_{t-2}}\right] $$ by fitting the model with $$ \hat{\gamma} = (R'(ZZ' \otimes \widehat{\Sigma_U}^{-1})R)^{-1}R'(Z' \otimes \widehat{\Sigma_U}^{-1})(\operatorname{vec}{Y} - (Z' \otimes I_k )r) $$ (equation 5.2.9), and computing the test statistic $$ \frac{(C_0\hat{\gamma}-c)'(C_0((ZZ')^{-1}\otimes\widehat{\Sigma_U})C_0')^{-1}(C_0\hat{\gamma}-c)}{N} \sim F_{N,\,T-Kp-1} $$ where $Z$ is the matrix of lags of $Y=(w, x, y, z)$, and $C_0$ is the matrix that restricts $\beta_A$ to $\beta_0$. And yes, I know this is equivalent to fitting many separate VARs, but I want to test for lag length simultaneously. I am using the Wald test rather than AIC, etc. at the suggestion of Toda & Yamamoto (1995), who find that its asymptotic distribution is a standard $\chi^2$ as long as the lag order is greater than or equal to the largest order of integration among all components of $Y$. And I'm assuming that none of these variables has more than one unit root, since it's a) implausible in this data and b) I don't have enough data points to distinguish a feature as complicated as an $I(2)$ structure from, say, a time-varying trend. **update:** I think I'm I just looking for a matrix operator that maps $\mathbb{R}^k$ to $\mathbb{R}^l$ where $l < k$. Is that accurate? Should I ask this on Math.SE?


title: Is there a value to show how fragmented is a group?

text: I would like to know if there is a way to calculate a number that can show how fragmented is a group? It should be something like "variance of categorical data"? For example I want to see how fragmented is Android/iOS in terms of OS version. The number would be high for Android, where there are a lot of versions, and the usage percentage of each version is small. For iOS, the number would be small, since there is a dominant usage percentage of one version.


title: glm with offset for rate data

text: I am attempting to fit a GLM to rate data. In my case it is the number of provisioning visits to a bird nest per hour. I model the data like in the example below including time as an offset. `number of visits ~ x , offset=log(time), family=poisson(link = log)` However this is complicated by the fact that my nests also have different numbers of chicks and I want to include this in my model. This suggests to me that I need to adjust my offset to account for both time and number of chicks. Am I thinking about this the correct way or do I need a different approach?


title: Given MTTF and Number of Items, how to calculate failing parts with Time?

text: I am wondering how to make use of MTTF Here is the situation, I am given an MTTF for an item type x and a certain demand for that item in the next 25 years, say 100 parts that will be in operation starting from now. I need to make sure that failing parts within the next 25 years get replaced, by storing more parts in inventory equal to the number of parts that will fail according to the MTTF. So given an MTTF of 10000 Hrs and 100 Parts, how many will fail within 25 years? Using a negative expon distribution?


title: Is Slice Sampling a special case of Gibbs Sampling?

text: I read on this thread the following: > If you _can_ use both the gibbs sampler and slice sampling to sample from a > posterior I would use the Gibbs sampler as the slice sampler seems > unnecessary to me <...> If you cannot use the gibbs but you can use the > slice then your question seems irrelevant. But then, when on these slides from M. Jordan's lecture he says that: > Slice sampling is a special case of Gibbs sampling ( ** _in a product > space_** ). What does he mean by that? Also, **when** would one not be able to do Gibbs sampling but still do slice sampling?


title: what is the meaning of the Samples in NER?

text: I would like to know in NER (Named Entity Recognition ) problem , which concept should be considered as samples? each token as a sample? or each sentence ? or each Named Entity should be considered as a sample?


title: How to plot data output of clustering?

text: I tried clustering a set of data (a set of marks) and got 2 clusters. I would like to graphically represent it. Bit confused about the representation, since I don't have the (x,y) coordinates. Also looking for MATLAB/Python function for doing so. **EDIT** I think posting data make the question clearer. I have two clusters I made using kmeans clustering in Python (not using scipy). They are               class 1: a=[3222403552.0, 3222493472.0, 3222491808.0, 3222489152.0, 3222413632.0,      3222394528.0, 3222414976.0, 3222522768.0, 3222403552.0, 3222498896.0, 3222541408.0,      3222403552.0, 3222402816.0, 3222588192.0, 3222403552.0, 3222410272.0, 3222394560.0,      3222402704.0, 3222298192.0, 3222409264.0, 3222414688.0, 3222522512.0, 3222404096.0,      3222486720.0, 3222403968.0, 3222486368.0, 3222376320.0, 3222522896.0, 3222403552.0,      3222374480.0, 3222491648.0, 3222543024.0, 3222376848.0, 3222403552.0, 3222591616.0,      3222376944.0, 3222325568.0, 3222488864.0, 3222548416.0, 3222424176.0, 3222415024.0,      3222403552.0, 3222407504.0, 3222489584.0, 3222407872.0, 3222402736.0, 3222402032.0,      3222410208.0, 3222414816.0, 3222523024.0, 3222552656.0, 3222487168.0, 3222403728.0,      3222319440.0, 3222375840.0, 3222325136.0, 3222311568.0, 3222491984.0, 3222542032.0,      3222539984.0, 3222522256.0, 3222588336.0, 3222316784.0, 3222488304.0, 3222351360.0,      3222545536.0, 3222323728.0, 3222413824.0, 3222415120.0, 3222403552.0, 3222514624.0,      3222408000.0, 3222413856.0, 3222408640.0, 3222377072.0, 3222324304.0, 3222524016.0,      3222324000.0, 3222489808.0, 3222403552.0, 3223571920.0, 3222522384.0, 3222319712.0,      3222374512.0, 3222375456.0, 3222489968.0, 3222492752.0, 3222413920.0, 3222394448.0,      3222403552.0, 3222403552.0, 3222540576.0, 3222407408.0, 3222415072.0, 3222388272.0,      3222549264.0, 3222325280.0, 3222548208.0, 3222298608.0, 3222413760.0, 3222409408.0,      3222542528.0, 3222473296.0, 3222428384.0, 3222413696.0, 3222486224.0, 3222361280.0,      3222522640.0, 3222492080.0, 3222472144.0, 3222376560.0, 3222378736.0, 3222364544.0,      3222407776.0, 3222359872.0, 3222492928.0, 3222440496.0, 3222499408.0, 3222450272.0,      3222351904.0, 3222352480.0, 3222413952.0, 3222556416.0, 3222410304.0, 3222399984.0,      3222494736.0, 3222388288.0, 3222403552.0, 3222323824.0, 3222523616.0, 3222394656.0,      3222404672.0, 3222405984.0, 3222490432.0, 3222407296.0, 3222394720.0, 3222596624.0,      3222597520.0, 3222598048.0, 3222403552.0, 3222403552.0, 3222403552.0, 3222324448.0,      3222408976.0, 3222448160.0, 3222366320.0, 3222489344.0, 3222403552.0, 3222494480.0,      3222382032.0, 3222450432.0, 3222352000.0, 3222352528.0, 3222414032.0, 3222728448.0,      3222299456.0, 3222400016.0, 3222495056.0, 3222388848.0, 3222403552.0, 3222487568.0,      3222523744.0, 3222394624.0, 3222408112.0, 3222406496.0, 3222405616.0, 3222592160.0,      3222549360.0, 3222438560.0, 3222597024.0, 3222597616.0, 3222598128.0, 3222403552.0,      3222403552.0, 3222403552.0, 3222499056.0, 3222408512.0, 3222402064.0, 3222368992.0,      3222511376.0, 3222414624.0, 3222554816.0, 3222494608.0, 3222449792.0, 3222351952.0,      3222352272.0, 3222394736.0, 3222311856.0, 3222414288.0, 3222402448.0, 3222401056.0,      3222413568.0, 3222298848.0, 3222297184.0, 3222488000.0, 3222490528.0, 3222394688.0,      3222408224.0, 3222406672.0, 3222404896.0, 3222443120.0, 3222403552.0, 3222596400.0,      3222597120.0, 3222597712.0, 3222400896.0, 3222403552.0, 3222403552.0, 3222403552.0,      3222299200.0, 3222321296.0, 3222364176.0, 3222602208.0, 3222513040.0, 3222414656.0,      3222564864.0, 3222407904.0, 3222449984.0, 3222352096.0, 3222352432.0, 3222452832.0,      3222368560.0, 3222414368.0, 3222399376.0, 3222298352.0, 3222573152.0, 3222438080.0,      3222409168.0, 3222523488.0, 3222394592.0, 3222405136.0, 3222490624.0, 3222406928.0,      3222407104.0, 3222442464.0, 3222403552.0, 3222596512.0, 3222597216.0, 3222597968.0,      3222438208.0, 3222403552.0, 3222403552.0, 3222403552.0]          class 2: b=[3498543128.0, 3498542920.0, 3498543252.0, 3498543752.0, 3498544872.0,      3498544528.0, 3498543024.0, 3498542548.0, 3498542232.0]      I would like to plot it. I tried the following and got the following result when I plot `a` and `b`.               pylab.plot(a,'x')     pylab.plot(b,'o')     pylab.show()      ![enter image description here](http://i.stack.imgur.com/CHIWH.png) **can I get a better visualization of clustering?**


title: Where can I obtain more weather data?

text: Where can I obtain more weather data? NOAA has some:   * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sflnd/DS.synop/   * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sflnd/DS.metar/   * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sfmar/DS.dbuoy/   * http://weather.noaa.gov/pub/SL.us008001/DF.an/DC.sfmar/DS.ships/ but I also want amateur observations. Specifically:   * Where can I get APRS WX data? I know that "telnet rotate.aprs.net 10152" with "user READONLY pass -1" will give me lots of information, but:      * Can I filter it down to just weather reports?      * Is this all APRS WX data, or just a subset?    * wunderground.com has over 17K stations at http://www.wunderground.com/weatherstation/index.asp reporting to them. Is this data available publicly on some "free weather exchange network"? I realize I can get it from wunderground.com, but would prefer a "legit" source.    * Are there other sources of weather data out there? Is there any one source that conglomerates multiple other sources? Recently found RAWS: http://raws.wrh.noaa.gov/rawsobs.html but sense there is more


title: Testing difference between two (adjusted) r^2

text: Let's say I have two regression models, one with three variables and one with four. Each spits out an adjusted r^2, which I can compare directly. Obviously, the model with the higher adjusted r^2 is the better fit, but is there way to test the difference between the two adjusted r^2 and get a p-value? I know you can do Chow test to test the difference between slopes, but this is variance, so I don't think that's what I'm looking for. Edit: One model does not simply contain a subset of variables from the other model, or else I would probably use stepwise regression. In model 1, I have four variables: W, X, Y, and Z. In model 2, I have three variables: W, X, and (Y+Z)/2. The idea is that if Y and Z are conceptually similar, the model may make better predictions by grouping these two variables together prior to entering them into the model.


title: How can I calculate the autocorrelation of a signal in Mathematica environment?

text: I tried `CorrelationFunction[Transpose[{data,data}]][[All,1,2]]` but it doesn't work! I mean the results are identical with those if I run `CorrelationFunction[data]`.


title: What are the disadvantages of the profile likelihood?

text: Consider a vector of parameters $(\theta_1, \theta_2)$, with $\theta_1$ the parameter of interest, and $\theta_2$ a nuisance parameter. If $L(\theta_1, \theta_2 ; x)$ is the likelihood constructed from the data $x$, the profile likelihood for $\theta_1$ is defined as $L_P(\theta_1 ; x) = L(\theta_1, \hat{\theta}_2(\theta_1) ; x)$ where $ \hat{\theta}_2(\theta_1)$ is the MLE of $\theta_2$ for a fixed value of $\theta_1$. $\bullet$ Maximising the profile likelihood with respect to $\theta_1$ leads to same estimate $\hat{\theta}_1$ as the one obtained by maximising the likelihood simultaneously with respect to $\theta_1$ and $\theta_2$. $\bullet$ I think the standard deviation of $\hat{\theta}_1$ may also be estimated from the second derivative of the profile likelihood. $\bullet$ The likelihood statistic for $H_0: \theta_1 = \theta_0$ can be written in terms of the profile likelihood: $LR = 2 \log( \tfrac{L_P(\hat{\theta}_1 ; x)}{L_P(\theta_0 ; x)})$. So, it seems that the profile likelihood can be used exactly as if it was a genuine likelihood. Is it really the case ? What are the main drawbacks of that approach ? And what about the 'rumor' that the estimator obtained from the profile likelihood is biased (edit: even asymptotically) ?


title: GLM model is singular

text: I am trying to fit a GLM model to my data (7 million rows, 153 variables) using R. More precisely I am using the Revoscaler package, but I suppose my issue would apply to other software as well. The calculation results in an error, saying that the GLM model is singular. I tried to find an explanation, without success. Could someone explain me what it means that the model is singular, and what needs to be done to solve the problem?


title: t-test vs. comparing regression coefficients

text: I understand that independent t-tests compare two means and establish whether they are statistically different from one another. If I want to compare whether there are statistical differences in the coefficients from two linear regressions, what kind of test would I use? Is this also a t-test?


title: Remove variance from temporal observations (input to HMM emission) by manipulating the individual observations

text: I have 52 weekly observations of count data that has some spikes and a definite trend (available here).               obs.pw <- read.csv(paste0(download.path,'obs-per-week.csv'))     obs.pw <- obs.pw[ 2:dim(obs.pw)[2]]      The mean weekly values calculated as such (can't post images yet):               plot(colMeans(obs.pw), xlab="Weeks", ylab="Mean Counts")      The weekly observations (8976 users per week) follow a zero-inflated negative binomial distribution. They are generally about 95% zeros. The summary data for the column sums of this data:               summary(colSums(obs.pw))     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      853    1424    1567    1607    1770    2374      Although the spikes and trend information is interesting it is not what I am interested in exploring now. I am **trying to the segment the non-zero counts for users into different ranges based on how infrequent that observation is**. I want to transform each observed integer count into ordinal categorical variables. Example variables would be: low, medium or high. A threshold is set for each category and the value is greater than or equal to the threshold the observation is assigned to that category (e.g. if(obs >= 30) cat.obs <- medium). The categories would ideally divide different quantiles of the data. The problem with the spread of points and trend highlighted in the graph means that **these thresholds are time dependent**. I am using these categorical variables as input to a Hidden Markov Model and want an observation to be just as rare on week 1 as week 50. I have tried dividing each individual observation by the mean observation count for that week but this line still seems to have a trend. Is there a way to scale these observations with the given mean data, or otherwise so that a threshold on day 1 for a low,med,high observation is roughly the same as for week 50? Or is this the wrong approach entirely?


title: Why doesn't this represent a normal approximation to the binomial?

text: Suppose the registrar's office at a college reports 58% of the students live on campus. An intern working in the administration building is unaware of this 58% parameter value. He designs a study in which he will take a random sample of 200 students and estimate the population proportion of all students that live in the college dorms using the resulting sample proportion. If he were to repeat his study many times, the possible values of the sample proportion would vary by about **__ away from the expected value of __** , on average. I understand the first blank is the standard deviation and the second is the mean, but why can't I approximate them using $\mathcal{N}(np, (np(1-p))^.5)$, which is the normal approximation to a binomial distribution? So that I have $\mathcal{N}(116, 6.98)$? So the first blank is 6.98 and the second is 116?


title: Finding the z score and p-value of a binomial distribution

text: Emily is a big fan of lady gaga, and 20% of the songs on her ipod are lady gaga songs. Suppose Emily has her ipod on shuffle and repeat mode, which can be assumed to mean that each song to be played always has the same 20% chance to be a lady gaga song. What is the approximate probability that at least 14 of 75 songs that Emily listens to will be from lady gaga? I understand I have to calculate the z score and find the corresponding p value with 1-P(z), but should I use the normal approximation to the binomial? Namely, $\mathcal{N}\left(np, \sqrt{np(1-p)}\right)$? I'm confused, because the answer lists $np$, and $np(1-p)$ instead of $\sqrt{np(1-p)}$. Where does this deviation come from?


title: Type of Distribution

text: Suppose a data set looks like exponential growth. See this. Would you use an exponential distribution to model this?


title: How can one speed up this correlation calculation in R without multicore?

text: I have a colleague who calculates correlations in which one set of scores for a subject (e.g. 100 scores) is correlated with another set of scores for that same subject. The resulting correlation reflects the degree to which those sets of scores are associated for that subject. He needs to do this for N subjects. Consider the following dataset:               ncol <- 100     nrow <- 100     x <- matrix(rnorm(ncol*nrow),nrow,ncol)     y <- matrix(rnorm(ncol*nrow),nrow,ncol)      The correct output vector of correlations would be:               diag(cor(t(x),t(y)))      Is there a faster way to do this without using a multicore package in R?


title: Confidence error bars and "central point": Should we emphasize the median?

text: Say I want to plot summary data with a point and a 95% confidence interval around that point. What should my point _really_ be? Mean, mode, or median? I know that mean = median for any symmetrical distribution, and the most common distributions used in data analysis (t distribution and normal distribution) have the convenient property of mean = mode = median, but what about the log-normal distribution: $\frac{1}{{x\sqrt {2\pi {\sigma ^2}} }}\exp \left( { - \frac{{{{\left( {\ln x - \mu } \right)}^2}}}{{{\sigma ^2}}}} \right)$ transformed back onto a linear scale? This is a common situation in biology (at least it would be if biologists were better with their statistics). Here are my arguments for all three moments. I cannot decide which is the strongest argument, so I am asking here:   1. **Median** The lower bound of the 95% confidence interval is the value which divides the lower 2.5% of the distribution from the upper 97.5%. Likewise, the upper bound divides the lower 97.5% from the upper 2.5%. Thus, the point that divides the upper 50% from the lower 50%, the median, should be the point between the upper and lower 95% CI bounds. In the case of the log-normal distribution, this is equal to $e^{\mu}$.   2. **Mode** If I am estimating a true value via regression, I want to show the value that is _the most likely_ of any value. That would be the peak of the distribution- the mode. In the case of the log-normal distribution, this is equal to $e^{\mu-\sigma^2}$, which is the lowest of the three ~~moments~~ measures of central tendency.   3. **Mean** Along the same logic used for mode, but now I am not just looking for the most right answer, but the answer that will minimize how wrong I'll be if I don't get it right. In the case of the log-normal distribution, this is equal to $e^{\mu+\frac{1}{2}\sigma^2}$, which is the highest of the three ~~moments~~ measures of central tendency.     * Which do you think is the strongest case?     * Is the answer the same if the value being plotted came from a regression model vs. summarizing raw data (note: I know that summaries of raw data are just single-parameter regressions, but not many biologists make that distinction).     * at the heart of the issue: what are we _really_ trying to show with the point in the center of a confidence limit? EDIT 01: @user603's answer has some good discussion and a request for more details. Here's some context that made me ask the question in the first place: Let's say I am doing regression on a dose-inhibition experiment with different drugs added to cultured cells. My model for the regression is: $ M = F+\frac{C-F}{1+\exp{\left( H(\ln{dose} - I_{\ln} \right) }}$ where $M$ is the measurement, $F$ is the lower asymptote (signal floor), $C$ is the uninhibited response (ceiling), and $I_{\ln}$ is the natural log of the half-inhibitory concentration for a particular drug. I do my regression under the assumption that $I_{\ln}$ is log-normal, so I get a regressed value (and confidence limits) for $I_{\ln}$. I want to report my findings in a nice chart showing the half-inhibitory concentration for the drugs on a linear scale. The goal of this or any scientific reporting of measure is to give our best guess at what is True in Reality, which we can approximate by experimentation and repeated measures. So I guess my question just became much more philosophical: what are summary statistics in science really trying to show? Our best estimate for a value (mode)? The estimate at which we have a 50% chance of over/undershooting (median)? The estimate by which we have the lowest probable deviation from the Truth (mean)? I often see it taught that confidence limits are are based on standard error (of the mean), but it is starting to seem to me that it should really be standard error (of the median) which just so happens to be (of the mean) when we're dealing with normal- and t-distributed uncertainty. So far, I am inclined to agree with @user603. Follow up question: If indeed I did want to show mean $\pm$ confidence limits, or mode $\pm$ confidence limits, how would those be calculated? Can you have standard error of the mode?


title: Statistical models cheat sheet

text: I was wondering if you know if a stastical model "cheat sheet(s)" that lists any or more information:   * When to use the model   * When not to use the model   * required and optional inputs   * expected outputs   * has the model been tested in different fields (policy, bio, engineering, manufacturing, etc)   * is it accepted in practice or research   * expected variation / accuracy / precision   * caveats   * scalability   * deprecated model, avoid or don't use   * etc .. I've seen hierarchies before on various websites, and some simplistic model cheat sheets in various textbooks; however, it'll be nice if there is a larger one that encompasses various types of models based on different types of analysis and theories.


title: How to evaluate quality of probability estimator for Bernoulli experiments?

text: Given that I have a set of bernoulli experiments, each with a different and unkown probability $p_i$ and an outcome $x_i$, and an estimator that for each experiment gives a prediction of the probability of the event, I want to measure the prediction quality of the estimator. Example: I have a stack of n "unfair" coins, each with a different probability $p_i$ for heads and $1-p_i$ for tails. The probabilities are unknown and I can flip each coin only once. Assume that there is a "coin flipping expert" which can have a close look at each coin before flipping them and make an estimate for the probabilities, based on form, size, width, regularity and so on. After the expert makes his prediction, the coin is flipped and the result is noted. After all coins are flipped, I want to measure how good the expert was, for example on a scale between 0 and 1, where 1 means perfect prediction and 0 means pure randomness. I would also be interested in bias / variance of the predictor.


title: Interpretation of logloss value

text: Does anyone have a interpretation of a logloss value? Am I correct to assume that values closer to 0 and 1 are more likely to be an indication that the predicted value is incorrect? Thanks.


title: Characterizing the inter-arrival time of software threads

text: For a multi-threaded application, I want to identify the nature of the application based on the arrival times of each thread. Example(Are thread launch spaced regularly, are they bursty in nature or are they highly dynamic). To identify this, the starting time of each thread is measured. The variance is calculated for time taken between every thread launch. This gives an idea of the spread in inter-arrival times. More spread indicates high dynamism. Based on the percentage of zscores that fall within the first level of standard deviation, a prediction is made on the level of dynamism (if more than 50% of values falls within one standard deviation - low dynamism otherwise high dynamism). Is this approach correct? Also, the values obtained for the starting time of each thread is only for a sample and not for the whole population. For this reason, the standard deviation is calculated using n-1. Are there any exisiting techniques for better estimation? Would clustering the values based on similarity help in any way? I'd highly appreciate it if you could point me to any reference. Thanks in advance.


title: Correlation with change score and baseline

text: My question is similar to this question, but I was not allowed to comment on it, so I will ask again tailored to my specific problem: I computed a correlation between a measure of brain volume and change in a standardized school achievement score (controlling for variables as age, IQ, gender) and get statistically significant results. However, I am not sure how to interpret these as I think that change in the achievement test is dependent on the baseline score (participants with better baseline scores are more likely to show decline, thoses with worse increases, i.e., regression to the mean). If I understand correctly, it is not appropriate to control for the baseline value. What else can I do to ensure that the found effects represent an association between structural brain difference and potential to change in achievement and are not just the result of regression to the mean? Would analyzing a subgroup without the participants who showed more extrem values at baseline be adequate?!


title: Are there any probabilistic models for graph-based recommender systems?

text: All I can find now is somehow based on random walks or graph kernels, which is nice, but I want to have a more or less solid probabilistic foundation for my recommender system for bounds and estimations and stuff which usually comes with probabilistic models. I probably should clarify: I have a graph and want to build recommender engine what recommends certain vertexes based on the ones selected by user.


title: Can the gamma distribution act as conjugate in this MVN case?

text: I have a sort of an expression for a multivariate normal-gamma distribution as follows: $$ P(\lambda|w) = P(w|\lambda) P(\lambda) $$ $P(w|\lambda)$ is a multivariate normal with 0 mean and parameterised with the precision matrix given by $(\lambda\Lambda)$. $\lambda$ is a scalar. $P(\lambda)$ is a gamma distribution with scale $s$ and rate $t$. I was under the impression that $\lambda$ would be gamma distributed. So, I have the following: $$ P(\lambda|w) \propto \lambda^{s-1} \exp(-t\lambda) \:|\lambda\Lambda|^{1/2} \exp(-\frac{w'(\lambda\Lambda)w}{2}) $$ Now, the problem is it does not seem possible to formulate this in the familiar Gamma form when one of the terms is a multivariate rather than a univariate Gaussian?


title: Central limit theorem versus law of large numbers

text: The central limit theorem states that the mean of i.i.d. variables, as $N$ goes to infinity, becomes normally distributed. This raises two questions:   1. **Can we deduce from this the law of large numbers?** If the law of large numbers says that the mean of a sample of a random variable's values equals the true mean $\mu$ as $N$ goes to infinity, then it seems even stronger to say that (as the central limit says) that the value becomes $\mathcal N(\mu, \sigma)$ where $\sigma$ is the standard deviation. Is it fair then to say that central limit implies the law of large numbers?   2. **Does the central limit theorem apply to linear combination of variables?** Thanks.


title: Seeking text describing the Mantel-Haenszel test

text: Seeking a textbook or other publication describing the Mantel-Haenszel test of a $2 \times 2 \times r$ table that provides at least one numeric example. Thanks.


title: Generating functions for Poisson regression using MATLAB

text: Alright, I'm having an insane amount of difficulty for what seems like a simple concept. I need to generate a bunch of basis functions for a curve that underlies some simulated Poisson distributed data and then run Newton's method on it to fit the log-likelihood (Poisson regression). I understand the math behind Newton's method enough to code that up and fit the data. However, I'm having a lot of difficulty just trying to generate a set of basis functions for which the equation $-\mathbf{X}^{\mathrm{T}}\cdot\mathrm{diag}\left(\exp\left(\mathbf{X}\cdot\boldsymbol\theta\right)\right)\cdot\mathbf{X}$ results in a tridiagonal matrix, where $$\mathbf{X}\cdot\boldsymbol\theta= \begin{bmatrix} f_{1}\left(\mathbf{t}_1\right) & f_{2}\left(\mathbf{t}_1\right) & \cdots & f_{d}\left(\mathbf{t}_1\right)\\ f_{1}\left(\mathbf{t}_2\right) & f_{2}\left(\mathbf{t}_2\right) & \cdots & f_{d}\left(\mathbf{t}_2\right) \\ \vdots & \vdots & \ddots & \vdots\\ f_{1}\left(\mathbf{t}_n\right) & f_{2}\left(\mathbf{t}_n\right) & \cdots & f_{d}\left(\mathbf{t}_n\right) \end{bmatrix} \begin{bmatrix} \theta_1\\ \theta_2\\ \vdots\\ \theta_d \end{bmatrix},$$ and $\mathrm{diag}$ results in a matrix with the elements of its argument on the diagonal of an otherwise zero- filled matrix. Here's what I've got so far:               function ddl = fastnewton(d, varargin)         %% boilerplate for optional function arguments         p = inputParser;              % required arguments         p.addRequired('d', @isposintscalar)              % optional arguments         p.addOptional('n', 1000, @isposintscalar)              % keyword arguments         p.addParamValue('order', 3, @isposintscalar)         p.addParamValue('tminmax', [-10, 10], ...                 @(x) (isvector(x) || iscell(x)) && length(x) == 2)         p.addParamValue('weightfunc', @(varargin) abs(randn(varargin{:})), ...                 @(x) isa(x, function_handle))              % parse the args and get the results         p.parse(d, varargin{:})              r = p.Results;         d = r.d;         n = r.n;         order = r.order;              if iscell(r.tminmax)             r.tminmax = cell2mat(r.tminmax);         end              tminmax = r.tminmax;         weightfunc = r.weightfunc;              %% make a time vector between tmax and tmin, inclusive         % tminmax == [-10, 10]         tmin = tminmax(1);         tmax = tminmax(2);         t = linspace(tmin, tmax, n); % n == 1000              %% make some random weights using the 'weightfunc' function         % weightfunc == @(varargin) abs(randn(varargin{:}))         theta = weightfunc(1, d);              %% make the bump funcs         breaks = linspace(tmin, tmax, d + order);         coeffs = theta;         pp = spmak(breaks, diag(coeffs));              %% compute the value of the function to be estimated         Fbig = spval(pp, t);         F = sum(Fbig);              %% initialize the sample vector         y = poissrnd(exp(F));              subplot(211)         hold on         plot(t, y, 'LineWidth', 1)         plot(t, Fbig, 'LineWidth', 2)         plot(t, F, 'r', 'LineWidth', 4)         hold off         legend({'y', '', 'sum(F)'})              %% plotting         lw = 3;         opts = {'LineWidth', lw};         subplot(212)         hold on         plot(t, y, 'b', opts{:})         plot(t, F, 'g', opts{:})         hold off         axis tight              X = Fbig';         ddl = sparse(-X' * diag(exp(F)) * X);              figure         hold on         plot(t, y, 'b', 'LineWidth', 1)         plot(t, F, 'r', 'LineWidth', 2);         hold off              figure         spy(ddl)     end      You can run this code as follows: `fastnewton(30, 'order', 2)`. That will give you 30 basis functions of order 2 which will give you a tridiagonal matrix. However, if you run `fastnewton(30, 'order', 3)` you get nice, smooth basis functions but a matrix with a bandwith of 3 instead, which is not tridiagonal. MATLAB's cute little `spy` function gives you a graphical representation of the sparsity of a matrix. All I really need is some way to adjust the basis functions so that they are smooth AND the equation above results in a tridiagonal matrix, then I can do the rest myself. Thanks!


title: Cross validation procedure - is this right?

text: Just want to check that I am performing my cross validation procedures right. I'm using a non-linear svm. I do a five fold cross validation (5 splits of test/train on my original training data) and for each fold, run a grid search to find the optimal parameters for that train/test pair (e.g. best parameters where model is fitted on the training data then evaluated on the test data). After five iterations, for each parameter set (2 hyper-parameters), i have 5 fit measures. I just use the average of those 5 #s and find the parameter set with the best average. Does this sound correct? I'm pretty new to this so am not entirely aware of what other methods are there. Additionally, was wondering a few more things: 1) Any way to make it faster? Cross validation + grid search is pretty computationally intensive. 2) Any other validation methods rather than k-fold (stratified I might add) cv? My data are timeseries so I was considering a moving window type validation, but wasn't sure. Any thoughts welcome. I should add that I'm asking because my training fit (the average over the 5 cv folds) is still much better than my actual test data fit. I'm trying to figure out what might be causing this and the best way to reduce this difference. I realize increasing the # of folds may help though it also raises question 1) as well. Thanks guys and thumbs up on the site. Great info.


title: How to compare two non-randomly selected groups within a sample?

text: I am quite new to Statistics, so please do not flame me for such a simple question. I made a controlled experiment and successfully gathered the data. 41 persons participated. They made a questionnaire and a game. From both of them, they obtained a single, discrete score. The situation is the following: I split the 41 people in two groups A and B (those that obtained a score <= 0 in the questionnaire vs. those that obtained a positive score in the questionnaire). Actually, A is composed by just 6 of them but different splits may be possible under different reasoning. Let's just keep this one for now. Therefore, I have two means for the game, X and Y, respectively. I would like to see if there are significant evidences that Y > X. Which test should I use? I am really confused. I continue to read around that if the sample size is > 30, I may use the T test anyway. Then, there is the Mann–Whitney U test to be used for two samples if we can not assume Normal Distribution but should be applied for two samples of independent observation. Any hints? Thank you very much. Update: Ideally speaking, the population could also be categorized by the questionnaire score, I am trying to see if there is statistical evidence that people with a positive score also would have a higher game score


title: Is there a version of the arcsine rule for biased coins?

text: I'm just reading about the arcsine rule for the first time in the context of coin flipping, i.e., that one player is very likely to lead a large percentage of the time in a finite number of trials. **Question** : Is whether there is a generalized version of this for biased coin? I tried finding something via a search engine, but a lot of the literature is in the context of academic papers on Brownian motion, so I'm unable to make heads or tails of it (pun unintended). Thanks!


title: How to cluster survey data?

text: I have designed a rather long (250 Qn) survey designed to uncover user clusters. The questions are such that the pattern of answering should elicit user clusters, but I am having trouble uncovering these with my analyses to date. For example, a typical Qn might be: 'Are you more of a dog or cat person?' Or, 'Chocolate or vanilla?' etc. I'm coding these questions in a binary format. So, if the user answered the two questions above with 1. Dog and 2. Vanilla, the user's answer matrix would look like: [1 0 0 1] Signifying that the user chose the first and fourth answer, where the answer space is [Dog Cat Chocolate Vanilla] I have roughly 300 respondents who have answered all 250 questions, giving around 800 possible answers, so my binary [user x answer] matrix is 300 x 800. I have run SVD on this matrix. The first factor relates to the number of people who selected that answer (magnitude) as expected. The second factor clusters nicely into male / female (I know because I ask gender) respondents. My problem is, all other factors are Gaussian and offer no way for me to split them into groups. A plot-matrix of the factors shows no grouping whatsoever. A clue: when I look at the highest and lowest factor values for factors 3, 4, & 5, I can determine that there are definite personality types represented. For example, cautious/risky or conservative/outgoing or frugal/outlandish. But these are just the tails of a Gaussian. I am completely unable to separate these questions by anything but a 'random' threshold to the Gaussian tail. The goal is to have a subset of the 250 questions that would allow me to quickly characterize a respondent, but right now, the only clustering I am able to assign is that of gender. Thanks in advance!


title: Comparing odds and payouts of different lotteries

text: I'm not a statistician, but I'm interested in how far I can push myself with my basic algebra skills, armed with a C# development background and Excel. Suppose I have several lotteries, each with different payouts, and odds, what is the right way to go about comparing and ranking the various games? What would a graph of that data look like? (Scatter with a logarithmic axis?)


title: What is an assignment matrix?

text: I'm trying to implement a topic model using a Latent Dirichlet allocation (LDA) algorithm. I'm using sentences as my dataset. What is Ck in the given instructions? The instructions are as follows: For dataset D: Run LDA with k topics on D to produce assignment matrix Ck Its on page 29 of this pdf http://sambrody.info/Lectures/Brody_aspect- sentiment_Google.pdf


title: overdispersion and zero

text: I have run GLM and GAM for my data (Mnemiopsis leyidi) in 4 season. I have used the quasipoisson for reduce of the effect of over-dispersion. in some the papers write that negative binomial is better than quasipoisson. Moreover, dependent variable in spring and winter have a number of zero. I really do not know that I should use which of the distributions. Can i use zero-infalted for quasipoisson? What method or index is used to determine the over-dispersion in GAM? if possible for you, please tell me appropriated distribution in each season for GLM in spring: n=132, number of zero=34, family=poisson Null deviance: 1779.13 on 131 degrees of freedom Residual deviance: 154.07 on 91 degrees of freedom AIC: 557.83 for glm in winter: n=132, number of zero=23, family=poisson Null deviance: 3063.05 on 131 degrees of freedom Residual deviance: 561.68 on 94 degrees of freedom AIC: 1056.8 for glm in summer: n=132, number of zero=11, family=poisson Null deviance: 7494.9 on 131 degrees of freedom Residual deviance: 3102.0 on 99 degrees of freedom AIC: 3780.2 for glm in autumn: n=132, number of zero=2, family=poisson Null deviance: 16016.5 on 131 degrees of freedom Residual deviance: 4922.9 on 96 degrees of freedom AIC: 5780.8


title: How many models I need?

text: I am doing an estimation using a bunch of sparse data. Suppose that I have a 100x100 grid and 20 data is available on this 2D grid. One solution is to use a determiastic method and estimate the other unknown locations on this grid, but as you may know, it does not guarantee the heterogeneity and we always can expect a very smooth results which is so far from the reality. Therefore, I am going to use the Monte-Carlo based estimation methods (such as sequential gaussian simulation) and produce a lot of estimation maps and then take an ensemble. My question is that how many model really I need to create? Because I can create an infinite number of models, but as it is clear, it is not affordable to do this.


title: Post hoc test for object of class glmerMod in R?

text: I am running generalised linear mixed effects models in R using the lme4 package. I am wondering if there are any post-hoc tests available for models built using the glmer function? I know the package lmerTest has a function for post hoc testing lmer models (class merMod), but it won't work for objects of class glmerMod (method given in the lmerTest manual). My models have either a numeric or integer response variable, two fixed effects (both a factor with 2 levels), and two random effects (both a factor, one nested within the other). I have used either a Gamma (for numeric response) or poisson (for integer response) error family. I want to know which combinations of levels of the fixed effects are statistically different to other combinations.   1. Are there any options for post hoc testing this type of model?   2. If not, can anyone recommend another method for performing generalised linear mixed effects models in R that do allow for post hoc testing?


title: Chi-squared calculation for (nonlinear) model goodness of fit for data with trials

text: I would like to be certain of how to calculate the $\chi^2$ value for this scenario: > We desire a measure of the goodness of fit of a nonlinear model to a data > set. The data consists of many data points $x_i$, where $i$ is an index, and > $x_i$ is specific to the experiment (e.g., measured in seconds, meters, > etc). At each data point $x_i$, a set of trials $\{y_{j}\}_x$ is collected. I frequently take data of this form. Usually, I simply average the trials and fit the results to an appropriate function. Until now, I've only needed to judge goodness of fit by eye, but now I need a more quantitative measure. Take for example the following data: ![enter image description here](http://i.stack.imgur.com/nRzID.png) There are 16 data points. At each point, 20 trials were "run", and the results are plotted above with black dots. The trials at each point are distributed normally, with a mean that follows a Voigt distribution across the entire data set, and with variances that are slightly different for each data point. The means are shown as red dots. What if we didn't know that the parent distribution of the data is a Voigt distribution? We might suspect that it is normally distributed. We could find the mean at each data point (red dots), and then fit them to a Gaussian and examine the goodness of fit. The goodness of fit may be measured by $\chi^2$, which is defined by Bevington & Robinson (2003) as $\chi^2 = \sum \frac{\left(y_i - y(x_i)\right)^2}{\sigma_i^2}$ where $y_i$ are the deviations of the data from the fitting function $y(x_i)$, and $\sigma_i^2$ are the variances at each point. My question is how to interpret this. Three different interpretations are: (a) At each data point, we sum the squared residuals $(y(x_i) - y_i)^2$, and divide the result by the variance for that point, $\sigma_i = \sum (y_j - \mu_i)^2 / (N_i - 1)$. This results in a list of $\chi^2$ values for each data point: ![enter image description here](http://i.stack.imgur.com/6V5Aq.png) I might then take the mean of these values and report the result as the $\chi^2$ statistic for the entire data set. For the above data, the resulting $\chi^2$ using this method would be 1.082. (b) I use the _mean_ at each data point for $y_i$ in the above definition. I again use the set of trials at each point to calculate a variance at each point (as in (a)). A pseudo-formula for $\chi^2$ would then be: $\chi^2 = \sum \frac{\left( \text{fit result at data point} - \text{mean of trials at data point} \right)^2}{\text{variance of trials at data point}}$ There is a question of whether the denominator should be the variance of the trials w.r.t the mean at that data point, or w.r.t the fit value. The following plot shows the residuals (red) and variances (gray). ![enter image description here](http://i.stack.imgur.com/t9Eqa.png) Assuming the former, the resulting $\chi^2$ for the above data is 1.243. (The latter gives 0.933.) (c) I consider the data as a set of scans where in each scan I have only one trial per data point. In other words, I have 20 independent sets of measurements of the data set distribution. (This is not unrealistic - data is often collected in this manner!) I fit each scan, and then calculate the $\chi^2$ value for each one, using the method of (b). The variance for each data point is calculated the same way as in (a) and (b): by considering all the trials at each data point. The result is 20 separate $\chi^2$ values, one for each "scan". I might then average these. Once $\chi^2$ is known, it must be interpreted. One possibility would be to convert it to the reduced chi-squared statistic by dividing by the degrees of freedom, and then comparing to one. This is problematic for nonlinear curves (see for example this paper). Another possibility is to use the probability distribution for $\chi^2$ to check if $P_{\chi}\left(\chi^2, \nu \right)$ is reasonably close to 0.5. Before I seriously consider either, though, I want to be sure that I'm calculating $\chi^2$ correctly.


title: Obtaining variance decomposition table with categorical variables in the model in R

text: I am running a multiple regression with the code:               model = lm(Price ~ factor(Make) + Cylinder + Liter + Doors + Cruise + Sound +leather, data=kuiper)      I got the `summary(model)` output without any error. I am able to get the output of variance decomposition table when I use the categorical variables in the model, but I am not able to produce the output for variance decomposition table (VD Table) using the `colldiag` function in package `perturb` when I have categorical variables in the model. R is throwing the following error:               print(colldiag(model))      Error in v^2 : non-numeric argument to binary operator      Is there any alternative R code for producing the variance decomposition table for categorical variables?


title: Can an estimator of the mean of a distribution with no variance have a variance?

text: Suppose you have a sample from a distribution with a mean but no defined variance, like the Pareto with tail parameter between 1 and 2, or Student’s t with 2 degrees of freedom. Can an unbiased estimator of the mean of such a distribution have a (defined, finite) variance?


title: Choosing a Sample Size for a SOM Cluster

text: Is there an accepted formula to determine a good lower bound for the number of samples needed for the initial input on an SOM clustering algorithm? For example - If I know that there are X number of attributes, and I suspect that there are at least Y distinct groups, and I have an input of 4 million line items - it will take a few days to process, so we'd like to take a simple random sample of those 4 million but we don't know how many we should take. A related question is, how can we predict how long it would take to process x number of entries with y number of unique attritubes.


title: Hypergeometric card question

text: I have a set of 60 things. Of those things, 24 belong to one type, 8 to another type, and 4 to a third type. If you select 7 things from the set, what's the probability of getting at least one of each together?


title: How to avoid overfitting when using crossvalidation within Genetic Algorithms

text: This is a long set-up, but the pure intellectual challenge will make it worthwhile I promise ;-) I have marketing data where there is a treatment and a control (i.e a customer gets no treatment). The event of interest (getting a loan) is relatively rare (<1%). My objective is to model the incremental lift between the response rate of the treatment and the control (treated group book rate - control group book rate) and use the model to make decisions about who to promote to in the future. The treated group is large (600,000 records) and the control is about 15% the size. This is a marketing exercise and we want to target those who need to be targeted to take the action of interest and not waste funds on those who will "do it anyway". I have hundreds of variables and have experimented with various forms of Uplift modeling AKA Net Lift Models. I have tried many of the state-of-the-art methods in the literature and common practice. None very stable on this data set unfortunately. I know (theoretically and after some experimentation) that there are a few variables that might impact the incremental lift. So, I created a matrix with the combinations of the levels of these variables and the number of records in the treated group, the number in the control group and the number of events of interest in each. So, from each row in the matrix one can calculate the incremental lift. There are 84 rows in the matrix. ![enter image description here](http://i.stack.imgur.com/As623.jpg) I was think of modeling this (difference in) proportion using a beta regression, but the counts in some rows are very spares (perhaps no records in the control and more frequently, there are no events of interest). This can be seen in the top couple rows of the sample data above. I began thinking about searching for the optimal solution to which of the rows of the matrix to select. Rows that are selected have the number of treatedHH and treatedLoans summed, along with the control. I am looking to maximize profit which can be estimated from these numbers. I pushed the data through a genetic algorithm to determine which rows to keep. I got a solution returned and the result was better than including everyone (which is the base case). But, when I ran that selection on the validation sample I partitioned, the result was not so. **My question** : Is there a way to design cross validation into this fitness function so that the solution does not over fit - which I presume happened in my first attempt. Here is the fitness function I used:               calcProfit<-function(selectVec=c())     {              TreatLoans<-sum(selectVec*dat$TreatedLoans)          ControLoans<-sum(selectVec*dat$ControlLoans)         TreatHH<-sum(selectVec*dat$treatedHH)          ControlHH<-sum(selectVec*dat$controlHH)                   Incre.RR<-(TreatLoans/TreatHH)-(ControLoans/ControlHH)         Incre.Loans<- Incre.RR * TreatHH         Incre.Rev <- Incre.Loans*1400         Incre.Profit<- (-1)*(Incre.Rev - (0.48*TreatHH))              Incre.Profit               }      and the call in R: `rbga.results = rbga.bin(size=84, zeroToOneRatio=3,evalFunc=calcProfit,iters=5000)`


title: How do I emphasize a speed increase between two sets of data?

text: I've run a test on a program, upgraded the program and then run the same tests again. The first test showed that the median time for the program to respond was approximately 2.0 seconds. The second test showed that, after the upgrade, the median time for the program to respond was approximately 0.5 seconds. While writing a paper, I really want to emphasize the _speed increase_ of the program after the upgrade. I recognize that I can multiply 0.5 seconds by 4 to get 2 seconds, but that really means, "the pre-upgraded program was 4 times slower than the upgraded program". I want to say something like, "After the upgrade, the program ran x times faster than before the upgrade." What is the proper way to get x?


title: R's equivalente of scikit's KFold

text: I'm new to R and I'm trying to set up a basic k folds CV loop. In Python I'd use scikit's KFold.               import numpy as np     from sklearn.cross_validation import KFold     Y = np.array([1, 1, 3, 4])          kf = KFold(len(Y), n_folds=2, indices=False)          for train, test in kf:           print("%s %s" % (train, test))     [False False  True  True] [ True  True False False]     [ True  True False False] [False False  True  True]      In R I was trying to use createFolds from the caret package, but it seems as though it's trying to do a stratified sample on the response variable.               library(caret)          Y <- c(1, 1, 3, 4)          fold <- createFolds(Y, 2)          fold     $Fold1     [1] 2 3          $Fold2     [1] 1 4      All I need is a simple k folds split. How do you do that in R?


title: Finding the bionomial probability?

text: Ok so i have been taught this formula regarding binomial probability > Repeat an even n times >   >   >       x = # of successes >       p = probability of successes >       q = probability of failure >   > > p(x = a) = nCa * p^a * q^(n-a) Ok wow seems good. So i try to solve this question using this formula. > A pair of fair dice is rolled 10 times. Let X be the number of rolls in > which we see at least one 2. > > What is the probability of seeing at least one 2 in any one roll of the pair > of dice? The probability of seeing at least one 2 = 1 - probability of seeing no 2s at all So i calculate the probability of seeing no 2s at all. x = 0 p = 1/6 q = 5/6 n = 10 Using the formula, p(x=0) = 10C0 * (1/6)^0 * (5/6)^10 Which gives me, p(x=0) = 0.1615 Probability of seeing at least one 2 = 1 - 0.1615 = 0.83 And i checked the ans, it stated the ans was 0.306. So where did i go wrong? If someone can help me please, i would be much grateful


title: Logistic regression for ranking: how do you represent inter-human variation?

text: This machine learning contest in Kaggle has a benchmark solution that uses logistic regression: http://www.kaggle.com/c/predict-who-is-more-influential- in-a-social-network The context provides data where two twitter users are compared based on the values of same features. So for each data row in the learning set we have               Twitter_User_A_More_influential | Twitter User_B_moreInfluential | {User_As_feature_vals} | {User_Bs_feature_vals}      If user A is more influential, the value is 1 (and B is 0) and the reverse if B is more influential. What is interesting is that the benchmark solution ( https://gist.github.com/fhuszar/5372873 ) is using logistic regression. I am not sure if logistic regression really encapsulates the model we have here, and I'd like to know how it may be improved (if it may be improved at all), especially to include inter-human variation. See the snippet from the python code used in the benchmark solution:               X_train = transform_features(X_train_A) - transform_features(X_train_B)     model = linear_model.LogisticRegression(fit_intercept=False)     model.fit(X_train,y_train)      This code transforms feature values of twitter user A and B (using new_ val = log(1+x)) and builds an array that has the difference between the transformed/regularized values. y_train contains ranking values for A: 1 or 0. The code fits logistic regression using training data and ranking values of A. The simplifications I see in this implementation are:   1. It turns ranking into classification, expressing more influential as influential or not.    2. It expresses the relation between the difference (change) in feature values and preferences of human agents.   3. It does not include anything that represents variation of judgement from N agents.  Number 3 feels like the biggest simplification to me. As it is, this model seems to represent a single individual's ranking of N comparisons, but the definition of data mentions more than 1 user (presumably close to N). I've been thinking about how to include that aspect in logistic regression, but so far I could not find a way of doing it. How would you represent the N human agents in terms of a random variable here? I'd appreciate different ways of modelling this data, from a statistics perspective.


title: How to reduce dimensionality of vector?

text: I have a set of vectors. I'm working on ways to reduce a n-dimensional vector to a unary value (one-dimensional), say (x1,x2,....,xn) ------> y This single value needs to be the characteristic value of the vector. Each unique vector produces a unique output value. Which of the following methods is appropriate:   1. norm of the vector - square root of sum of squares that measures euclidian distance from origin   2. compute hash of F, using some hashing techniques avoiding collision   3. use linear regression to compute, y = w1*x1 + w2*x2 + ... + wn*xn. Unlikely to be good if there is no good dependence of input values on output   4. feature extraction technique like PCA that assigns weights to each of x1,x2,..xn based on the set of input vectors.


title: Performing model validation in stata

text: I need to validate a cox-regression model internally - i was thinking of either k-fold cross validation or bootstrapping methods. are either of these possible to do with a cox model? any idea how to perform these? an example of our stepwise regression is below: sw, pr(0.05): stcox age3 apricat2 meldcat2 etiol5a sex if yrs>=0.5 Thanks for your help


title: Correlation between numbers versus percentages

text: Suppose we have three variables $X_1, \dots, X_3$. There are $5$ categories for each variable. In each of the categories, the variables can take at most $3$ values. Suppose we are interested in the frequency of values in each category across the $3$ variables. Suppose the data is the following: $X_1 = (1,2,3,3,2)$, $X_2 = (2,3,1,3,1)$ and $X_3 = (1,3,1,3,1)$. These represent the observed counts. To determine the correlations between each variable should we look at the counts? Or should we look at percentages? In other words, should we look at $X_1 = (1/3, 2/3, 1,1, 2/3) \dots$?


title: Comparing similarity of two fingerprint or bitvector

text: I have a question that's two part. 1) The first part is that I have a bit vector or a fingerprint of a digital circuit where the each bit is high if a specific structure or component exists and zero other wise. All the fingerprint is of equal length. The bit position indicates a specific structure and is the same throughout all of them. So for example, bit 1 = 2bit mux, bit 2 is a 3 bit mux, bit 4 is a 4 bit mux and so on. Using something like Tanimoto compares it well, but it is a bit too strict. For example, if I have b1 and b2 where                              X    V     b1 = 0000 1000 1111 1000     b2 = 0001 1001 0011 0010                  X        V      we can see where there are exact matches(where the 1's line up), however, the bit position where there is an X or V over it i also want to be able to take that into account as well because even though they are not over each other, their bit position is close enough that it is "similar" to it. Are there any methods or algorithms that could possibly take this into account? Would a sequence alignment be the right path to take, or should it be more of a distance metric using something like a sliding window? 2)The second part is similar to the fingerprint, however, instead of bitvectors, it is an ordered set where each position represents how many of a specific structure exists. So in this case, I have               s1 = 157 523  708 0 0 19 0 0  0  0  0     s2 = 148 537  711 0 0 15 0 0  0  0  0      s3 = 32  139  124 0 0 15 0 0 16 16 16      Is there a way to tell how similar s1 is to s2 and s3? I would expect s1 to be similar to s2 and not s3. The order of the numbers are important since it indicates how many of a specific component exists. I have looked into cosine similarity but the results are not especially matching what i expect. Thank you


title: Compare linear regression slopes for two datasets, each with a different predictor variable

text: I have two data sets. In one of them I am controlling the predictor variable $X_1$ and observe the response on $Y$. In the second one I control $X_2$ and observe the response on $Y$. I perform linear regression for each data set separately. $X_1$ and $X_2$ are different by 6 orders of magnitude. In theory, are the estimated slopes between the two data sets comparable? I suppose the F-test corrects for the difference in orders of magnitude? (There may be a correlation between $X_1$ and $X_2$, but I think in my case it is okay to assume their correlation is zero.)


title: Doing low-dimensional KNN on a large dataset

text: I have a simple two-dimensional dataset with columns X1,X2, and [outcome], and I want to try KNN (probably K around 100 or 1000, though ideally CV would be possible). The problem is that my dataset has a couple million rows. Are there any preferred packages/approaches for dealing with this sort of thing?


title: Help with R and grouping/aggregate/*apply/data.table

text: I'm very new to R and having trouble running functions to get the answers I need. I have example data PCSTest http://pastebin.com/z9Ti3nHB that looks something like this:               Date        Site            Word     --------------------------------------     9/1/2012    slashdot        javascript     9/1/2012    stackexchange   R     9/1/2012    reddit          R     9/1/2012    slashdot        javascript     9/1/2012    stackexchange   javascript     9/5/2012    reddit          R     9/8/2012    slashdot        javascript     9/8/2012    stackexchange   R     9/8/2012    reddit          R     9/8/2012    slashdot        javascript     9/18/2012   stackexchange   R     9/18/2012   reddit          R     9/18/2012   slashdot        javascript     9/18/2012   stackexchange   R     9/27/2012   reddit          R     9/27/2012   slashdot        R      My goal is to looks for trends in the occurrences of different words as they relate to sites over time. I can count them:               library(plyr)        PCSTest <- read.csv(file="c:/PCS/PCS Data - Test.csv", header=TRUE)     PCSTest$Date <- as.Date(PCSTest$Date, "%m/%d/%Y")     PCSTest$Date <- as.POSIXct(PCSTest$Date)     countTest <- count(PCSTest, c("Date", "Site", "Word"))      which gives this:                                 Date          Site       Word freq     1  2012-08-31 20:00:00        reddit          R    4     2  2012-08-31 20:00:00      slashdot javascript    7     3  2012-08-31 20:00:00 stackexchange javascript    1     4  2012-08-31 20:00:00 stackexchange          R    2     5  2012-09-01 20:00:00        reddit javascript    2     6  2012-09-01 20:00:00      slashdot          R    3     7  2012-09-04 20:00:00        reddit          R    1     8  2012-09-07 20:00:00        reddit          R    1     9  2012-09-07 20:00:00      slashdot javascript    2     10 2012-09-07 20:00:00 stackexchange          R    1     11 2012-09-09 20:00:00 stackexchange javascript    4     12 2012-09-10 20:00:00      slashdot          R    4     13 2012-09-14 20:00:00        reddit javascript    4     14 2012-09-17 20:00:00        reddit          R    4     15 2012-09-17 20:00:00      slashdot javascript    1     16 2012-09-17 20:00:00 stackexchange          R    2     17 2012-09-19 20:00:00        reddit javascript    2     18 2012-09-23 20:00:00 stackexchange javascript    2     19 2012-09-24 20:00:00        reddit javascript    3     20 2012-09-24 20:00:00 stackexchange javascript    1     21 2012-09-24 20:00:00 stackexchange          R    4     22 2012-09-25 20:00:00        reddit javascript    5     23 2012-09-25 20:00:00      slashdot javascript    3     24 2012-09-25 20:00:00 stackexchange          R    7     25 2012-09-26 20:00:00        reddit          R    1     26 2012-09-26 20:00:00      slashdot          R    5      or plot them all:               library(ggplot2)     ggplot(data=countTest, aes(x=Date, y=freq, group=interaction(Site, Word), colour=interaction(Site, Word), shape=Site)) + geom_line() + geom_point()      ![My plot of Frequency per day for Words per Site](http://i.stack.imgur.com/epGfu.jpg) I need to do some calculations on the data now, so I've tried aggregate               aggregate(freq ~ Site + Word, data = countTest,  function(freq) cbind(mean(freq), max(freq)))[order(-agg$freq[,3]),]      which gives:                          Site       Word freq.1 freq.2     2      slashdot javascript   3.25   7.00     5      slashdot          R   4.00   5.00     1        reddit javascript   3.20   5.00     4        reddit          R   2.20   4.00     6 stackexchange          R   3.20   7.00     3 stackexchange javascript   2.00   4.00      What I would like in that last result is a column that has the average frequency per day, something like... sum(freq)/20 days, calculated from the data, maybe a moving average even. Also, I would like another column with the slope/linear regression. How would I calculate that in the aggregate function? Or, how would I make any of this better/faster? I know that there are apply and data.table functions but I don't see how I would use them. Any help would be greatly appreciated!


title: Estimates of regression coefficients at subject- vs group-level

text: For each of several subjects, I have mean reaction times for each of several stimuli. For each individual subject, I want to regress the RTs (as Y) against the stimuli (as X) and see if I get a slope that is significantly different from zero. For that, I used Matlab's fitlm function which gives me, for the X predictor, a coefficient estimate and a p-value. I then want to see if, at the group level, if those coefficients (slopes) are significantly different from zero, using a 1-sample t-test. My questions: 1) For each subject, is it correct to look at the p-value as an indication of whether the slope is significantly different from zero, or should I look at a goodness-of-fit measure instead? I thought that, as in the case of correlation, the question of slope is more to do with the notion of effect size than that of statistical significance 2) Is it correct to say that those p-values represent within-subject variability while the p-value of the 1-sample t-test represents between- subjects variability? In other words, can the effect (significant linear dependence between Y and X) be present at the subject level but not necessarily at the group level? Thanks!


title: Predict observation using Hidden Markov Models

text: I have a sequence of observations e.g. ["Click","Scroll","Hover","Zoom","Select"]. I need to predict the next value of this observation sequence but not the next hidden state. I know that there are three fundamental problems for HMMs: a)Given the model parameters and observed data, we can estimate the optimal sequence of hidden states. b)Given the model parameters and observed data, we can calculate the likelihood of the data. c)Given just the observed data, we can estimate the model parameters. So, for solutions to my kind of problem:   * I thought by referencing to b) that if I make conditional sequences of data each of them ending with one of the possible values that could stand for the next observation and calculating the likelihood of each of them given the model, then can this be considered as prediction?  To be more specific, in my example (if I know that the possible observations can only be Click/Scroll/Hover/Zoom/Select) I will simulate the following sequences ["Click","Scroll","Hover","Zoom","Select","Scroll"] ["Click","Scroll","Hover","Zoom","Select","Click"] ["Click","Scroll","Hover","Zoom","Select","Select"] ["Click","Scroll","Hover","Zoom","Select","Hover"] ["Click","Scroll","Hover","Zoom","Select","Zoom"] and the sequence that gives higher probability is "the predicted", so eventually I also have the predicted next observation, which will be the last observation of the sequence that gives higher likelihod. Is this correct?   * Another way as it is referred in this link would be to predict the most likely hidden-state-sequence based on a) and then through the emission distribution of the last hidden state to calculate the mean of this distribution? The above link was never verified and I am wondering if anyone could verify it.   * Other way would be to get the sum of the likelihoods of all states each of them multiplied with the mean of the state's distribution. Is it correct? Thank you in advance for any feedback you can give me.


title: How to find a mapping function from n dimensional space to m dimensional space

text: What are different methods/approaches to find a function $T:x_1 \to x_2$ that minimizes a cost function defined as: $\text{cost}(T) = | f(Q(x_1)) - f(P(x_2)) | = | f(Q(x_1)) - f(P(T(x_1)) |$ where, $x_1$ has $N$ dimensions, $x_2$ has $M$ dimensions, and $P$ and $Q$ are any function from $M$ and $N$ dimensional space to $Z$-Dimensional space, respectively. $f$ is a function from $Z$-Dimensional space to 1-dimensional space. $N$, $M$, and $Z$ can be any positive integer.


title: Estimating training values in machine learning

text: I've started to learn machine learning using the book 'Machine Learning' by Tom Mitchell. The first chapter introduces the reader to machine learning and to the first problem: learning to play checkers. It starts by choosing a target function $V: Board \rightarrow \mathbb{R}$, which maps a particular board setting to a score (how good the current board state is). A score of $100$ means the game was won, $-100$ it was lost, $0$ draw. It then chooses a representation for the board: features. $$V = w_0 + x_1*w_1 + x_2*w_2 + ...$$ where $w_i$ are weights and $x_i$ are features (for example the number of pieces you currently have on the board). Features can be determined from the board configuration itself, but the weights have to be adjusted by the training examples. Each training example is a pair of the form $\langle b, V_{train}(b)\rangle$, where $b$ is the board state described by its features and $V_{train}(b)$ is the training value for $b$. Determining $V_{train}(b)$ for boards that are won or lost are easy, but we also need training examples for intermediate board states. This is where I start being confused. It states a good approximation for intermediate boards is $V_{train}(b) \leftarrow V(successor(b))$, where $successor(b)$ denotes the next board state following $b$ for which it is again the program's turn to move (i.e., the board state following the program's move and the opponent's response). How am I supposed to know what's $successor(b)$? Which move am I supposed to make and how do I know what will be the opponent's move?


title: GEV of Normal Distribution and relationship of the parameters

text: My question goes on Extreme Value Theory for the Normal distribution (www.math.ethz.ch/~embrecht/RM/chap7.pdf):   1. Which type of GEV (Generalized Extreme Value) distribution does the Normal distribution belong to (Weibull/Gumbel/Frechet)?   2. If we have the Normal parameters given ($\mu, \sigma$), can we infer the parameter values of its GEV?


title: How do decision tree learning algorithms deal with missing values (under the hood)

text: What are the methods that decision tree learning algorithms use to deal with missing values. Do they simply full the slot in using a value called missing? Thanks.


title: Final Model Prediction using K-Fold Cross-Validation and Machine Learning Methods

text: Similar threads: Feature selection for "final" model when performing cross-validation in machine learning Choosing a predictive model after k-fold cross-validation * * * My question is quite simple and is definitely related to the similar threads above but what I am looking for is a concrete yes/no to the question below: I am working on a regression problem where I have a target function of 1 variable that I am trying to predict using 5 explanatory variables. I have 1200 examples of the response and explanatory data. I decide to split my 1200 examples into a calibration set of 1000 examples and a test set of 200 examples. The calibration set is used to train my model and the test set is completely independent. Let's say I am using a Neural Network of a particular configuration/parametrization and I am looking to find the best possible network weights and biases such that it provides the best performance on my test set. To do this I have chosen to perform k-fold cross-validation on the calibration data. Let's say I opt for 10 folds. I thus produce 10 different calibrated models (using the training and validation sets for each k-fold) each of the same configuration using the Neural Network described above. I now want to use the Neural Network to provide an output on my test set using the parameters (weight and biases) determined from the k-fold cross-validation. **To produce the estimates on the test set do I simply average the weights and biases from each of the 10 different calibrated models and use this parametrization to produce outputs to compare with my test set for the target function?** Thank you everyone for their help!


title: Forecasting Hourly Time Series based on previous weeks and same period in previous year/s

text: **The Problem** I have been tasked with a similar problem to that described in Forecasting hourly time series with daily, weekly & annual periodicity. My data shows the number of times that one of 384 different events has occurred per hour between 01/01/2005 and 31/12/2013. I have been asked to provide a 28 and 84 day forecast showing the number of each type of event expected per hour. **My Experience** Before I go on I will clarify that I am a bit of a newbie in R (and statistics) and while I realise that this is not at all a simple task (nor possibly one that R is best suited to) I don't have access to any other tools and there is no one else in my organisation that is prepared to tackle this problem. **What I have tried so far** So I have restricted myself to just one event type and have tried using `auto.arima` and `HoltWinters` on the entire 75,000 hourly datapoints using timeseries `ts(data,frequency=8760)`, `ts(data,frequency=8760/12)`, and `ts(data,frequency=8760/52.1xxxx)`. I excluded `ets` because it doesn't cope with frequencies above 24. The results from these tests were pretty much flat line which isn't suitable for what I want to use the forecast for. Because each week is closely related to the previous week/s I tried to simplify things by only forecasting based on the previous 4 and 8 weeks using `ts(data, frequency=168)` (treating it as weekly data). The results from the `HoltWinters` experiment look promising and I have discarded ARIMA because it is still producing a straight line. Over a 28 day period in 2013 the total of the `HoltWinters` "forecasted" values based on the previous 4 weeks is within 10% of the actual values which is close but I think I could get closer. (This is probably not the best way to test the values but it seems like a good indicator.) **My Planned next steps** I am now wanting to try and factor in the same period for the year/s previous. That is, if I am forecasting for December 2013 I will use November 2013 and December 2012, 2011, etc. **The Questions:**   1. How do I forecast using the previous weeks and the same period for the previous years?   2. Is there a better way of doing this? (I am open to suggestions here, I am really only getting by, by testing and refining my method and am conscious of the fact that there is probably much better ways of doing this.)   3. It is quite likely that occurrences of different event types are linked (i.e. if event _a_ is happening a lot event _b_ will likely be happening a lot too, and event _c_ will probably not happen). As a next step is it possible be find and then factor in these correlations? _Sorry that I have not been able to attach my data but the computer system at work does not allow us to access places like dropbox and I cant see any other way of attaching a csv to this question. If it is required to answer my questions then I will take a copy home and upload it from there._


title: Observational study with many data Vs. Experimental study with fewer data: Which has more value?

text: So I have recently been debating two possibilities for determining the effectiveness of, let us say for simplicity, a single marking strategy for getting return customers. Let us assume that we have years of historical data where a company has, at certain uncontrolled times and at uncontrolled store locations, marketed to in-store customers in an attempt to get them to return (let's say they give out pamphlets at the check-out register, for example). My question is: What's a more effective strategy.   1. Using the historical data from times in the past when the store has marketed to try to estimate the effectiveness of the marketing to get customers to return. (e.g. marketed customers are 50% more likely to return in a week)   2. Using a designed experiment to give some randomly selected customers pamphlets and do not give them to the rest (control). Perhaps do this for a period of time (maybe a week) and track them for a while (maybe a month) to estimate the effectiveness of the marketing to get customers to return. The argument for option 1 is that you have much more data; however, you have confounding variables and no obvious control group to compare what would have happened if no marketing was done. The argument for option 2 is that you essentially have a Randomized Controlled Experiment; however, you have less data which takes time (money) to gather in a controlled fashion. I understand there may not be a concrete answer, but is there any theory or practical statistical reasoning that suggests when one or the other option is better. Just as a note: assume customers can be tracked for return visits and try to avoid arguments based solely on logistical issues (unless you feel they are very important).


title: What is an indicator in statistics and why is it used to refine models?

text: In a problem we are asked to "refine the fitted model by using an indicator for the outlier". What does it mean to use an indicator for an outlier?


title: Temporal Correlation in R

text: Having used stack overflow for years without signing up I now truly need some help! My data consists of samples collected from the same location in three different streams between 2012-2014. The data collected includes continuous phosphorus, nitrate and suspended sediment measures. Rainfall data from the area has also been downloaded. The samples for each parameter have been taken daily although there are holes in the data set so I am guessing I cannot call these regularly spaced data points. I was simply modelling using lme on log transformed data, although this was causing NaNs. I then realised I needed to incorporate temporal correlation at each site due to the daily auto-correlation and not just produce a nested design. Having spend two days working through error msgs resulting from using templates in "Mixed Effects Models and Extensions in Ecology with R", I am still unsure which design to use for the data I've collected. Advice would be welcomed on the design my model should be taking....thank you!


title: Error when using glmulti with aicc?

text: Does anyone know the solution to the error I get when I run glmulti with the aicc criteria: Error in if (length(lesCrit) == confsetsize && minouN - minou >= -deltaM && : missing value where TRUE/FALSE needed It goes away when the confsetsize is low (~5). Thanks!


title: Using non-significant data for modelling/predictions

text: **Two medical interventions A and B were compared within randomised controlled trial:**   * differences in progression free survival (PFS, i.e. how long patient lived with controlled disease, from entering study to progression) significant with p<0,001,   * differences in overall survival (OS, i.e how long patient lived, from entering study to death) non-significant with p=0,12. **We would like to model/predict PFS and OS beyond time of RCT. Which model should be chosen:**   1. joint modelling of OS (we don't allow differences between A and B, single OS function for two interventions), separate modelling of PFS (two different functions for A and B),   2. separate modelling of OS and PFS (four different functions for A and B, two for PFS, two for OS). Models are build by analysing goodness-of-fit (AIC/BIC, visual inspection and comparison with external data sources) to Kaplan-Meier curves. We consider using Weibull, log-normal, Gompertz, log-logistic and Gamma distributions. * * * Please note that English is not my first language.


title: Problem with parallel analysis with psych

text: I have a data set with several hundred variables and some thousand records. I'm reviewing the different ways for running a Principal Component Analysis and choosing the principal components. First I used the function `prcomp()` to get the advantage of using SVD. The first component explains more than 50% of the variance and the second a % more. The remaining components are very small. This is the figure I got: ![prcomp](http://i.stack.imgur.com/rw2fl.png) Then, by reading some posts here in Cross Validated I found that a recommended method for choosing the principal components was the so called "parallel analysis" and 2 libraries were recommended including the psych. So I decided to try that method. This is the code I used:               > fa.parallel(my.dataframe,               fa="PC",                n.iter=100,               how.legend=FALSE,                main="My Plot")      I received the following results in the command prompt:               Loading required package: MASS     In smc, the correlation matrix was not invertible, smc's returned as 1s     In smc, the correlation matrix was not invertible, smc's returned as 1s     The determinant of the smoothed correlation was zero.     This means the objective function is not defined.     Chi square is based upon observed residuals.     The determinant of the smoothed correlation was zero.     This means the objective function is not defined for the null model either.     The Chi square is thus based upon observed correlations.     In factor.stats, the correlation matrix is singular, an approximation is used     In factor.scores, the correlation matrix is singular, an approximation is used     I was unable to calculate the factor score weights, factor loadings used instead     Parallel analysis suggests that the number of          factors =  54  and the number of components =  44      Warning messages:     1: In cor.smooth(R) : Matrix was not positive definite, smoothing was done     2: In cor.smooth(R) : Matrix was not positive definite, smoothing was done     3: In cor.smooth(r) : Matrix was not positive definite, smoothing was done     4: In factor.stats(r, loadings, Phi, n.obs = n.obs, np.obs = np.obs,  :       In factor.stats, the correlation matrix is singular, and we could not calculate           the beta weights for factor score estimates     5: In cor.smooth(r) : Matrix was not positive definite, smoothing was done      And this is the resulting figure: ![fa.parallel](http://i.stack.imgur.com/8t8Iy.png) None of the PCA figures I had reviewed look like this, and I'm totally lost at analyzing its significance. I guess I should choose the first 44 principal components. I'd really appreciate if somebody could explain this to me a little bit.


title: Correlations of correlations using 3 data sets

text: I'm interested in how 14 environmental parameters co-vary with other environmental parameters and whether they show the same patterns in 3 regions. So 1st, I made a pairwise correlation table of the 14 environmental variables for each region. Then, I made another pairwise correlation table of correlation values from 3 regions (i.e. region1 vs region2 & region1 vs region3 & region2 vs region3). I found a good correlation (a similar pattern) between region2 and region3 (or called "sector" in graphs). ![Correlations](http://i.stack.imgur.com/uEwXv.png) What I want to do from here is to put the all correlations values from the 3 regions together in one graph and also calculate correlation values (i.e. correlations of correlation values) because I want to know which environmental variables are correlated well in all 3 regions. I thought of using PCA, but it's not really what I imagined to be the best way to present a pattern. Maybe a 3D graph would be good for presentation. I'm a `R` user. If you could instruct me with `R` language on how to make that kind of graph and also how to calculate 3-way correlations, that would be much appreciated. Thanks!


title: Test Cox proportional hazard assumption (Bad Schoenfeld residuals)

text: Using R I generated a Cox model looking like this               > summary(cox5 <- coxph(surv_ob_2 ~ age + AMP + nation + GESL + year, data=sd)) # now gender makes sense     Call:     coxph(formula = surv_ob_2 ~ age + AMP + nation + GESL + year,          data = sd)            n= 35800, number of events= 31873                                        coef  exp(coef)   se(coef)       z             Pr(>|z|)         age                    -0.0136536  0.9864392  0.0005184 -26.336 < 0.0000000000000002 ***     AMPofficial            -0.6058868  0.5455904  0.0338201 -17.915 < 0.0000000000000002 ***     AMPmarginal employment  0.5032381  1.6540687  0.0155578  32.346 < 0.0000000000000002 ***     AMPfreelancer          -0.7277312  0.4830036  0.0342688 -21.236 < 0.0000000000000002 ***     nationFormerYugoslavia  0.2758670  1.3176726  0.0304818   9.050 < 0.0000000000000002 ***     nationGermany           0.2354548  1.2654842  0.0341056   6.904     0.00000000000507 ***     nationother             0.1421815  1.1527859  0.0144441   9.844 < 0.0000000000000002 ***     GESLM                   0.1089281  1.1150822  0.0114915   9.479 < 0.0000000000000002 ***     year                    0.0526784  1.0540906  0.0009662  54.521 < 0.0000000000000002 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                                 exp(coef) exp(-coef) lower .95 upper .95     age                       0.9864     1.0137    0.9854    0.9874     AMPofficial               0.5456     1.8329    0.5106    0.5830     AMPmarginal employment    1.6541     0.6046    1.6044    1.7053     AMPfreelancer             0.4830     2.0704    0.4516    0.5166     nationFormerYugoslavia    1.3177     0.7589    1.2413    1.3988     nationGermany             1.2655     0.7902    1.1837    1.3530     nationother               1.1528     0.8675    1.1206    1.1859     GESLM                     1.1151     0.8968    1.0902    1.1405     year                      1.0541     0.9487    1.0521    1.0561          Concordance= 0.623  (se = 0.002 )     Rsquare= 0.17   (max possible= 1 )     Likelihood ratio test= 6690  on 9 df,   p=0     Wald test            = 6141  on 9 df,   p=0     Score (logrank) test = 6285  on 9 df,   p=0      I now want to test if the proportional hazard assumption is met. Therefore I run               > print(test.cox5 <- cox.zph(cox5, transform=rank)) ## Gender does not meet ph-assumption                                rho   chisq             p     age                     0.1099  435.85 0.00000000000     AMPofficial             0.1001  312.68 0.00000000000     AMPmarginal employment  0.0169    9.26 0.00233833597     AMPfreelancer           0.0332   34.78 0.00000000369     nationFormerYugoslavia  0.0171    9.31 0.00227422735     nationGermany           0.0105    3.50 0.06133812325     nationother            -0.0226   16.14 0.00005894624     GESLM                   0.0157    7.88 0.00500199726     year                    0.1047  226.84 0.00000000000     GLOBAL                      NA 1193.78 0.00000000000      These results look very bad to me. Do I have to assume that only nationGermany (hardly) meets the ph-assumption, or is there another explanation as well? My teacher supposed that a possible reason could be the large sample size ( _n=35800_ ). Is he right? What else can I do to check if the ph-assumption is met and to improve the model? If you know any literature covering this issue, I would be interested as well. **[edit]**   The data I analyse comes from a register based labor market database. I try to model the length of employments. The database covers several decades, therefore also macroeconomic circumstances could have changed with time. This is why I added age _and_ year in the model. I am however unsure if this is a good idea...


title: Optimize number of layers and neurons with an optimization algorithm

text: I have a neural network that i want optimize number of hidden layers and neurons in every layer using an optimization algorithm like `Imperialist competitive algorithm`. As you know i should set limits of every variable for optimization. I can have maximum 2 layers. Besides this my system select best composition of features (among 21 features) based on cost function of my neural network (5-fold cross validation on test data) using that optimization algorithm. So my optimization algorithm has `21+2` variables. `21` for features and `2` for number of neurons in every layer. What is your idea and what do you recommend for boundary of number of neurons in every layer? I'm using `patternnet` in `MATLAB`. As you know we have difference composition of features and neural structure in calling cost function by optimization algorithm.


title: Multinomial logistic regression and marginal effects

text: I am trying to calculate the marginal effects of a multinomial logistic regression. To do this I use the `mlogit` package and the `effects()` function. Here is how the procedure works (source : `effects()` function of `mlogit` package) :               data("Fishing", package = "mlogit")     Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")     m <- mlogit(mode ~ price | income | catch, data = Fish)     # compute a data.frame containing the mean value of the covariates in the sample     z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean),       catch = tapply(catch, index(m)$alt, mean),      income = mean(income)))     # compute the marginal effects (the second one is an elasticity     effects(m, covariate = "income", data = z)     effects(m, covariate = "price", type = "rr", data = z)     effects(m, covariate = "catch", type = "ar", data = z)      I have no problem with first step (`mlogit.data()` function). I think my problem is in the specification of the multinomial regression. My regression (for example with three variables) is on the form: `Y ~ 0 | X1 + X2 + X3`. When I try to estimate the marginal effects for a model with 2 variables, there is no problem, however for 3 variables R console returns me the following error: "Error in if (rhs% in% c (1, 3)) {: argument is of length zero " (translation from error in R console in french). To understand what is my problem I tried to perform a multinomial regression of similar shape on the dataset "Fishing", i.e.,: `mode ~ 0 | income + price + catch` (even if this form has no "economic" sense.) Again the R console returns me the same error for 3 variables but manages to estimate these effects for a model with two variables. This leads me to think that my problem really comes from the specification of my multinomial regression. Do you know how I could find a solution to my problem? Or could you suggest another logit multinomial regression form ? Thank you for your help :)


title: Comparing group of people/countries on multiple binary variables

text: I have a large amount of data from three countries on work values. People are given 17 work value items (good pay, working with pleasant people, job security etc.) and asked the ones that they find important. The ones that they find important are coded with 1 and the ones that are not mentioned re coded with 0. Now, I want to classify these work values into three types: intrinsic (8 work values), extrinsic (6 ww), social (4ww). Then I would like to see how they look in each country and most importantly, I would like to make comparisons among the countries. I am thinking of computing the items and doing a chi-square testing, but I don't feel confident about this solution. What do you think? Am I totally wrong? Do you have any other suggestions? I would appreciate if you could share your views on this. Thank you, Mel


title: Orthogonal intersection of linear family and exponential family

text: I asked the following question in MSE for which I couldn't get any answer yet. I thought this would be a better place for that question. In statistical maniolds $S=\{p_\theta\}$,$\theta=(\theta_1,\dots,\theta_n)$, the Riemaanian metric usually defined is the Fisher information metric $$g_{ij}(\partial_i,\partial_j)=\int \partial_i(\log p_\theta) \partial_i(\log p_\theta)~p_\theta~dx$$ The associated connection coefficients are defined by $$\Gamma_{ij}^k=\int \partial_i\partial_j(\log p_\theta)\partial_k (\log p_\theta)~ p_{\theta}~dx$$ where $\partial_i=\frac{\partial}{\partial\theta_i}$. My question is, what is the intuition behind defining these? Is there a way to **prove using the above metric and connection** that the linear family of probability distributions $$L=\{p:\int f_i(x)p(x)~dx=m_i, i=1,\dots,k\}$$ intersects "orthogonally" the associated exponential family $$\mathcal{E}=\{p:p(x)=c(\theta)q(x)\exp(-\sum_{i=1}^k\theta_i f_i(x))\}$$ in the sense that $L\cap\mathcal{E}=\{p^*\}$ where $p^*$ satisfies $$D(p\|q)=D(p\|p^*)+D(p^*\|q)$$ for every $p\in L, q\in \mathcal{E}$. I recently came to know about the connection between Fisher information metric and the relative entropy: $$D( p(\cdot , a+da) \| p(.,a) )\approx\frac{1}{2} g_{i,j} da^{i} da^{j}$$ Would this be a backbone in establishing the above result?


title: Analysing overdispersed data with generalised linear models

text: Let's say I have an explanatory variable and a response variable that represents counts. I want to see if the explanatory variable can predicts counts. I'm aware the response variable is overdispersed. What I probably should do is analyse the data with a generalised linear model using the negative binomial distribution. But let's say I ignore the overdispersion and analyse the data using poisson regression. What would be consequences of analysing these data with a poisson distribution? Would consequences be the same as ignoring heteroscedasticity in linear regression therefore incorrect standard errors/P-values?


title: Comparing two time series over different periods of time

text: I have two time series that I wish to compare their movements. The issue is that one time series is daily data since 1920 and the other is daily data since 2008. When plotting them separately over their time domains, I can see that they follow a very similar shape. My aim is to transform one of the time series to be on the same time scale as the other to determine whether they correlate or not. What is the best way to go about this? Thank you in advance for any suggestions!


title: Conditional distributions of ARCH/GARCH Model

text: I am using ARIMA-GARCH/ARCH modeling on financial time series. There are few querries regaring the general concepts. Like What is meant by a conditional distribution? Is it a distribution of stock returns (original data) or the distribution of the residuals which we get from ARIMA (conditional mean model)?


title: How can I make 2-way ANOVA (Type I, II and III) with R?

text: I have the following data with factor 1 (A, B and C) and factor 2 (D and E): $$ \begin{array}{ccc} \hline & D & E\\ \hline A & 68 & 59\\ & 65 & 57\\ & 63 & 54\\ & 59 & 56\\ & 67 & \\ \hline B & 59 & 51\\ & 50 & 45\\ & 51 & 46\\ & & 48\\ & & 49\\ \hline C & 40 & 47\\ & 39 & 39\\ & 35 & 40\\ & 36 & 40 \end{array} $$ Now I would like to make two-way ANOVAs of type I, II and III with R and interpret them. I made the following: ` rm(list=ls(all=T)) data <\- c(68,65,63,59,67,59,50,51,40,39,35,36,59,57,54,56,51,45,46,48,49,47,39,40,40) f1 <\- c(0,0,0,0,0,1,1,1,2,2,2,2,0,0,0,0,1,1,1,1,1,2,2,2,2) ### A=0,B=1,C=2 f2 <\- c(3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4) ### D=3,E=4 summary(aov(data ~ f1*f2)) ### type 1 install.packages("car") library(car) Anova(lm(data~f1*f2),type=c("II")) ### type 2 Anova(lm(data~f1*f2),type=c("III")) ### type 3 ` * * * For type 1 I get ` Df Sum Sq Mean Sq F value Pr(>F)   f1 1 1941.9 1941.9 193.420 4.58e-12 *** f2 1 65.5 65.5 6.527 0.01845 *   f1:f2 1 145.9 145.9 14.537 0.00102 ** Residuals 21 210.8 10.0   Signif. codes: 0 ‘ _**’ 0.001 ‘*_ ’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ` * * * For type 2 I get ` Response: data Sum Sq Df F value Pr(>F)   f1 1901.11 1 189.355 5.607e-12 *** f2 65.53 1 6.527 0.018450 *   f1:f2 145.95 1 14.537 0.001016 ** Residuals 210.84 21   Signif. codes: 0 ‘ _**’ 0.001 ‘*_ ’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ` * * * For type 3 I finally get ` Response: data Sum Sq Df F value Pr(>F)   (Intercept) 1802.06 1 179.489 9.308e-12 *** f1 329.87 1 32.856 1.089e-05 _**f2 208.54 1 20.771 0.0001714 **_f1:f2 145.95 1 14.537 0.0010157 ** Residuals 210.84 21   Signif. codes: 0 ‘ _**’ 0.001 ‘*_ ’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ` Is that allright? How would you interpret these results?