title: Obtaining $\lambda$ for a Poisson Distribution

text: Is it possible to manipulate a Poisson Distribution's $\lambda$ value such that it has a probability mass function that looks like the curve below with data points $(8,0)$,$(7,0.01)$,$(5,0.1)$ and $(2,1)$. I have been trying to obtain a $\lambda$ that will satisfy these values but have been struggling as the probability mass function equation does not break down easily. ![enter image description here](http://i.stack.imgur.com/z57Tj.jpg)


title: Getting a posterior belief from noisy observations

text: Say I have a dataset $X = \{x_i\} \supset \mathbb{R}^n$ for which I assume that $x_i = y_i + \sigma(y_i)$ for some unobserved variables $\{ y_i \}$. That is, I believe that my data is subject to heteroscedastic noise. From what I see, this rules out methods such as factor analysis. I am currently searching for methods for the estimation of $p(y_i|x_i)$ on my data. Does this situation have a name? Are there any canonical approaches?


title: Estimate the variance of MLE

text: From Wikipedia > Under the Wald statistical test, the maximum likelihood estimate > $\hat\theta$ of the parameter(s) of interest $\theta$ is compared with the > proposed value $\theta_0$, ... In the univariate case, the Wald statistic is > $$ \frac{ ( \widehat{ \theta}-\theta_0 )^2 }{\operatorname{var}(\hat \theta > )} $$ ... where $\operatorname{var}(\widehat\theta)$ is the variance of the > maximum likelihood estimate. My question is about how to estimate $\operatorname{var}(\widehat\theta)$, i.e. the variance of the MLE $\hat{\theta}$?   1. The same article says > A reasonable estimate of the variance for the MLE can be given by > $\frac{1}{I_n(\widehat\theta)}$ , where $I_n$ is the Fisher information of > the parameter.   2. I am thinking of another way. Since $\operatorname{var}(\widehat\theta)$ is $f(n, \theta, \eta)$, i.e. some function of the sample size $n$ and the true value $\theta$ of the parameter of interest and some other parameter $\eta$ (because $\theta$ alone may not be fully able to specify the distribution of the sample), can we first estimate the MLEs of $\theta$ and $\eta$ as $\hat{\theta}$ and $\hat{\eta}$, and then will $f(n, \hat{\theta}, \hat{\eta})$ be the MLE of $\operatorname{var}(\widehat\theta)$ because of the equivariance of MLE?   3. Will the above two methods be the same (in some cases, such as when the statistical model of the sample is normal distributions?) Thanks and regards!


title: Classification of signals, and Wavelet Transform based denoising

text: Stationary Wavelet Transform based denoising methods are considered the best because of their translational invariance. On the other hand, Matlab promises SWT based denoising but in the ddencmp function for denoising a signal, it is actually using Discrete Wavelet Transform when it should rather threshold the SWT co-efficients. Is there something that I am missing? I will be absolutely grateful.


title: Determining the "dominant" color in an image

text: I'm attempting to make a simple mosaic generator, where it takes a photo, then uses hundreds of other photos to represent each pixel. To do so, I want to determine the dominant color of each picture I am using as a tile, so I can know what colors it can represent. I have an array containing the RGB values of each pixel in the image, so it is not hard for me to manipulate the data. The issue is that the images I will be using are completely random, so I have no prior knowledge of the distributions of colors (meaning that simply using mean/mode etc would be very inaccurate). I am thinking some sort of specialized mode, since it might capture concentrated colors well, but I am not too sure. Any help would be appreciated. Thanks


title: How to find the connection between variables

text: Assume that you want to approximate your joint distribution on set of variables via decomposing it into several smaller distributions. Assume that you don't know any prior knowledge about the statistical or algebraic connection between variables; the only input is a set of samples from variables. How do you claim whether some variables tend to be dependent/independent? (note that this is more than just calculating correlation, which is indeed a linear operation.)


title: Attempting to investigate a relationship between one variable and another that has three possible values

text: I'm looking to investigate the relationship between a series of variables on a recurring basis. One variable will always be the dependent variable, and will take three possible values: `good`, `bad`, or `other`, which I have coded as `1`, `-1`, and `0`, respectively. I'm generally interested in seeing how a series of additional independent variables interact with this first variable, all of which will range linearly. For instance, one such variable might be called `minutes` and will range from zero to a few thousand. Essentially, I want to see how variables like `minutes` might affect the likelihood of a `good` outcome. Ordinarily, I would run a regression, but the fact that the dependent variable here takes three possible forms is making that tricky. What's the correct way to run this sort of a series of tests?


title: Calculating probability based on mixed variables

text: Assume there are $K$ people and iid. parameters $a_1,\ldots,a_K$ associated to them with $a_i \sim U(0,1)$. Person $i$ observes his own fixed $a_i$ with some noise: \begin{equation} X^{(1)}_i= a_i+ e^{(1)}_i, \end{equation} where the $e_i$ is an error terms distributed normally with $N(0, v)$. He learns $X^{(1)}_i$. After observing $X^{(1)}_i$, person $i$ would like to calculate the probability that next time the persons $1,\ldots,K$ get signals (denote the new signals by $X^{(2)}_i$), his signal will be greater than everyone else's signals, conditional on the value $X^{(1)}_i$. Everyone's signal has an error term with the same distribution, and error terms are independent. In summary, I would like to calculate the following: \begin{equation} \textrm{Prob}(X^{(2)}_i > X^{(2)}_{j}, \textrm{for all } j\neq i\mid X^{(1)}_i). \end{equation} If someone can help me with how to approach this problem, I would really appreciate it.


title: equating problem with instrument that have two version

text: I want to measure metacognition (as a whole construct and per subscale, (knowledge of cognition (KC)) and regulation of cognition (RC)) development in children using the Jr MAI inventory. The Jr MAI has two version, version A (for 9-11) and B (for 13-15), version A has 3 point Likert scales and version B has 5 point Likert scales and they also have different numbers of questions: Version A has 6 items of KC /RC while B has 8 items of KC / RC. Therefore how would I compare the mean (of group A (10-11) and B (13-14)) to see the developmental changes (metacognition/KC/RC), is there any standard method with reference?. My second question: Can different types of scores of different variable (e.g. raw score, composite likert scale score, standard score from standard manual test) be used in subsequent analysis such as regression, correlation or even in multiple mediator analysis (dependent and mediators)?.


title: Differentiating experimental data

text: I'm trying to differentiate experimental data which are measured with time. However, the time increments are not equally spaced. I first obtained a trend line for these data and then did a linear interpolation to obtain equally spaced points. Then I used central differencing to obtain the derivative. However, this gives a derivative that is not smooth and has steps in it. How can I solve this problem?


title: Logistic regression algorithm in Ruby

text: I have been using R to calculate logistic regression with many independent variables for a Ruby on Rails web application. However, I can no longer import data from the database to R using RPostgreSQL. The web host has stopped allowing insecure connections to the database. The point is, I either need to get a new web host, or write my own logistic regression algorithm in Ruby. Ruby probably isn't the best programming language for that kind of thing, but I don't really have a choice. Is there an easy to implement algorithm for multiple logistic regression?


title: What's the logic/motivation behind pre-filtering data based on variance?

text: I am currently working with gene expression data where I have a number of genes (variables) measured over a number of samples. I could not make sense of the statistics and visualization suite we have, so I started checking out the tutorial and came up with the following bit: > You specify a lower cutoff for the variance. You can either enter a value, > or use the slider. The value is relative to the variance of the variable > with the largest variance. The variance is calculated over the active > samples. > > Example: If you select the value 0.1, then all variables that have a > variance σ greater than or equal to 0.1 ·σ_max are kept, where σ_max is the > variance of the variable with the largest variance. I realize that it might be very fundamental question, but why would I want to filter my variables based on the variation alone? What I would want to know is which variables "explain" the difference between the samples (which are _hypothesized_ to belong to a couple of different biological conditions)


title: Significance of the slope of a straight line fit

text: I know the significance of the slope of a LMS linear regression can be calculated using the r2 coefficient of determination and looking up the appropriate value in an F table. However, I was thinking of perhaps making this "more robust" by replacing the LMS linear regression with a repeated median straight line fit, and perhaps even replacing the average value used to calculate r2 with the median value of the data. Is there any reason why this would NOT be a valid approach? Maybe the values in an F table are predicated on using LMS and averages, for example?


title: Next steps after obtaining a significant Pearson coefficient

text: Lets say I have an R value of .5 and therefore an $R^2$ value of .25, and a p value < 0.001. (All through the stock cor.test function in R). Is there a way to isolate those cases/observations whose variation is most explained by my $R^2$? My goal is two fold, find the group for which this correlation is maximized, and take the other group too explore separately for other interesting interactions instead if this one. Sorry for the imprecise wording here, clearly I'm not quite an expert in this area, and hence dont quite know how or what question to ask, but thank you for any direction you can provide.


title: A non-technical visualization of my database (a layman's ER Diagram?)

text: I would like to provide an intuitive visualization of my database for readers and users of the database that does not require technical knowledge of database structure. The figure will be included in a journal article, and the audience is primarily students and scientists, many of whom will not understand concepts such as "many to many" or keys. My goal is to convey the contents of the database with some insight into these relationships, but to convey a conceptual understanding rather than a technical one. A technical and comprehensive description of the database will be provided as an appendix, in the documentation. I see that there are many types of ER diagram, but I am not sure which one, if any is best suited for explaining a data model to my audience. How can I represent this in a non-technical way that avoids keys and provides a more intuitive understanding of the relationships among tables? Here is an example of my starting point: ![enter image description here](http://i.stack.imgur.com/RXsTU.png) The actual database will have 11 tables, including 3 used to define many-to- many relationships, and many fields that can be dropped.


title: Need help with normalizing data in SPSS

text: I apologize in advance... there's a good amount of stuff to read before you get to my question. Maybe I could have saved a word or two here and there... but overall, the background (I think) was required in order to best illustrate the topic. Ok, I recently purchased SPSS in support of an upcoming research analysis. I've not used SPSS before and would like to get some guidance/pointers for normalization. BACKGROUND: As part of this academic research, I plan to release a survey in the next few weeks. For my dependent variable (and dependent variable only), I will first provide the survey respondents a 2-column matrix where the first column lists the "name of a process" (across n number of rows) and the second column includes a checkbox. Let's say I have 10 default rows (i.e., processes) from which the participants check any processes (in column 2) in which they believe their team contributes to. Example: Survey participant #1 may check rows/processes: 1, 3, 7 Survey participant #2 may check rows/processes: 1, 4 Survey participant #3 may check rows/processes: 5, 7, 8, 9, 10 and so forth... Using branching logic in the electronic survey, their selected processes (e.g., 1, 3, 7) will be passed into subsequent questions. Their themes are as follows: Question #1 will require to indicate a frequency for how often a process may have changed. A survey participant will have to select a value -- on a scale of 1 through 7 -- for each process. Question #2, #3, and #4 will ask something else (note: for the purpose of this blog question, it's not important what the actual questions read). For the 2nd through 4th question, survey participants will have to provide answers of either "Yes" or "No". (... where 'yes' will be coded as "1"; and 'no' will be coded as "0" in SPSS). So far so good? ** Break ** Upon completion of the survey (I'm striving for several hundred participants), the data will be uploaded in SPSS. Now, again, all questions above are pertaining to a single dependent variable. Based on the survey participants' answers, I'd like to create the following: a) Individual score b) Overall organizational score Obviously, multiplication across the columns will not work... as soon as, e.g., as a participant selected "no" in either the 2nd, 3rd, and/or 4th question, that process would be scored as zero. I don't want that! So, addition is probably the better choice. Example of how survey participant #1 may have scored: 1: | 5 | 1 | 0 | 1 | = 7   3: | 2 | 1 | 1 | 0 | = 4   7: | 1 | 0 | 0 | 1 | = 2   Thus, participant #1's score = 13 (7 + 4 + 2) On the other hand, survey participant #3 may have scored: 5: | 3 | 1 | 0 | 0 | = 4   7: | 4 | 0 | 1 | 0 | = 5   8: | 2 | 0 | 0 | 0 | = 2   9: | 7 | 1 | 1 | 1 | = 10   10: | 3 | 1 | 0 | 1 | = 5 Thus, participant #3's score = 26 (4, 5, 2, 10, 5) Now, participant's #3 example score of 26 is twice as high as #1's example score of 13. This is primarily due to the larger number of projects/processes that the third survey participant is involved in. As I'm trying to measure some form of "health" of the organization, a larger total score (derived through a larger number of processes) could be viewed as a contributing factor to the "organization's health". In this example, this would be wrong! In other words, only larger scores on the individual process-level [(1 vs. 3 vs. 7) or (5, 7, 8, 9, 10) may be indicative of an "unhealthy" process/organization. So, process #9 by itself (with a score of 10) might be a contributor. So, my question is as follows: In SPSS, how can I normalize the data so that it won't make a difference as to how many processes a survey participant has selected (and thereby skewing the data)? Should I merely divide by the number of processes? So, 13/3 = 4.333 vs. 26/5 = 5.2 If so, how can I accomplish that in SPSS? I hope the above makes sense. Thank you, EEH


title: Difference between Taxonomy and Classification with respect to data/text mining

text: I am trying to understand the exact difference between Taxonomy and Classification with respect to data/text mining. I would appreciate it if you could explain with a simple example.


title: How to compute HDR probability level for Highest density regions estimates

text: In the `R` package `hdrcde` by Rob Hyndman for Computation of highest density regions in one and two dimensions, the hdrbw function Estimates the optimal bandwidth for 1-dimensional highest density regions. In the given example, HDRlevelVal variable is a probability 0 HDRlevel : HDR-level as defined in Hyndman (1996). Setting ‘HDRlevel’ equal to p (0 I have not found the referenced article. Can anybody give me some information about that article and function. Is there in `R`, packages that implements this function?               HDRlevelVal <- 0.55     x <- faithful$eruptions     hHDR <- hdrbw(x,HDRlevelVal)     HDRhat <- hdr.den(x,prob=100*(1-HDRlevelVal),h=hHDR)     Rhat <- hdr.den(x,prob=100*(1-HDRlevelVal),h=hHDR)


title: How do I calculate the range that 90% of the values fall into

text: Given a series of values, I know that 68% of the values fall withing one standard deviation and that 95% fall withing 2 standard deviations, but how do I calculate the range whereby 90% percent of the values will fall. Lets call that 90% range the typical range.


title: How to choose initial values for nonlinear least squares fit

text: The question above says it all. Basically my question is for a generic fitting function (could be arbitrarily complicated) which will be nonlinear in the parameters I am trying to estimate, how does one choose the initial values to initialize the fit? I am trying to do nonlinear least squares. Is there any strategy or method? Has this been studied? Any references? Anything besides ad hoc guessing? Specifically, right now one of the fitting forms I am working with is a Gaussian plus linear form with five parameters I am trying to estimate, like $$y=A e^{-\left(\frac{x-B}{C}\right)^2}+Dx+E$$ where $x = \log_{10}$(abscissa data) and $y = \log_{10}$(ordinate data) meaning that in log-log space my data looks like a straight line plus a bump which I am approximating by a Gaussian. I have no theory, nothing to guide me about how to initialize the nonlinear fit except perhaps graphing and eyeballing like the slope of the line and what the center/width of the bump is. But I have more than hundred of these fits to do so instead of graphing and guessing, I would prefer some approach which can be automated. I can't find any references, in the library or online. The only thing I can think of is to just randomly choose initial values. MATLAB offers to choose values randomly from [0,1] uniformly distributed. So with each data set, I run the randomly initialized fit a thousand times and then pick the one with the highest $r^2$? Any other (better) ideas? * * * Addendum #1 First, here are some visual representations of the data sets just to show you guys what kind of data am I talking about. I am posting both the data in its original form without any kind of transformation and then its visual representation in log-log space as it clarifies some of the data's features while distorting others. I am posting a sample of both good and bad data. ![good data](http://i.stack.imgur.com/xNCmg.jpg) ![log-log good data](http://i.stack.imgur.com/9EAA6.jpg) ![bad data](http://i.stack.imgur.com/4vAXC.jpg) ![log-log bad data](http://i.stack.imgur.com/iDJ34.jpg) Each of the six panels in each figure show four data sets plotted together red, green, blue, and cyan and each data set has exactly 20 data points. I am trying to fit each one of them with a straight line plus a gaussian because of the bumps seen in the data. The first figure is some of the good data. The second figure is the log-log plot of the same good data from figure one. The third figure is some of the bad data. The fourth figure is the log-log plot of figure three. There is much more data, these are just two subsets. Most of the data (about 3/4) is good, similar to the good data I showed here. Now some comments, please bear with me as this might get long but I think all of this detail is necessary. I'll try to be as concise as possible. I had originally expected a simple-power law (meaning straight line in log-log space). When I plotted everything in log-log space, I saw the unexpected bump at around 4.8 mHz. The bump was thoroughly investigated and was discovered in others work as well so its not that we messed up. It is physically there and other published works mention this too. So then I just added a gaussian term to my linear form. Note that this fit was to be done in log-log space (hence my two questions including this one). Now, after reading the answer by Stumpy Joe Pete to another question of mine (not related to these data at all) and reading this and this and references therein (stuff by Clauset), I realize that I shouldn't fit in log-log space. So now I want to do everything in pre-transformed space. ~~Question 1: Looking at the good data, I still think that a linear plus a gaussian in pre-transformed space is still a good form. I would love to hear from others who have more data-expreience what they think. Is gaussian+linear reasonable? Should I only do a gaussian? Or an entirely different form?~~ ~~Questions 2: Whatever the answer to question 1, I would still need (most likely) non-linear least squares fit so still need help with the initialization.~~ The data where we see two sets, we very heavily prefer to capture the first bump at around 4-5 mHz. So I don't want to add more gaussian terms and our gaussian term should be centered on the first bump which is almost always the bigger bump. We want "more accuracy" between 0.8mHz and around 5mHz. We don't care too much for the higher frequencies but don't want to ignore them completely either. So maybe some sort of weighing? Or B can be initialized around 4.8mHz always? The abscissa data is frequency in units of millihertz, denote it by $f$. The ordinate data is a coefficient we are computing, denote it by $L$. So no log transformation, and the form is $$L=A e^{-\left(\frac{f-B}{C}\right)^2}+Df+E.$$   * $f$ is frequency, is always positive.   * $L$ is a positive coefficient. So we are working in the first quadrant.   * $A$, the amplitude should always be positive too I think because we are just dealing with bumps. When I look at the data I always see peaks and no valleys. It looks like in all of the data there are multiple bumps at higher frequencies. The first bump is always much larger than others. In good data, the secondary bumps are very weak but in bad data (panel 2 and 5 for example), the secondary bumps are strong. So we actually _don't_ have a valley, but rather two bumps. Meaning that the amplitude $A>0$. And since we care mostly about the first peak, all the more reason for $A$ to be positive.   * $B$ is the center of the bump and we always want it on that big bump around 4-5mHz. In our resolved frequency range, it almost always appears at 4.8mHz.   * $C$ is the width of the bump. I imagine it is symmetric around zero meaning $C$ would have the same effect as $-C$ because it is squared. So we don't care what its value is. Let's say we prefer it positive.   * $D$ is the slope of the line, seems like it could be slightly negative so not enforcing any restriction on it. The slope is interesting in its own right so instead of enforcing any restrictions on it, we just want to see what it would be. Is it positive or negative? How big/small is it magnitude-wise? and so on.   * $E$ is the (almost) $L$-intercept. The subtle thing here is that because of the gaussian term, $E$ is not quite the $L$-intercept. The actual intercept (if we were to extrapolate to $f=0$) would be  $$Ae^{-(B/C)^2}+E.$$ So the only restriction here is that intercept should be positive as well. The intercept being zero, I don't know what that would mean. But negative would seem nonsensical for sure. I guess here we can allow for $E$ to be slightly negative with a small magnitude if necessary. The reason $E$ and the intercept is important here but some of our colleagues are actually interested in extrapolation as well. The minimum frequency we have is 0.8mHz and they want to extrapolate between 0 and 0.8mHz. My naive idea was to just use the fit to go all the way down to $f=0$. I know extrapolation is harder/more dangerous than interpolation but using a straight line plus a gaussian (hoping it decays fast enough) seems reasonable to me. Kind of like natural cubic splines with natural boundary conditions, the slope at the left endpoint, just extend the line and see where it crosses the $L$ axis. If it isn't negative then use that line for extrapolation. ~~Questions 3: What do you guys think extrapolating this way in this case? Any pros/cons? Any other ideas for extrapolation? Again we only care about the lower frequencies so extrapolating between 0 and 1mHz...sometimes very very small frequencies, close to zero. I know this post is already packed. I asked this question here because the answers might be related but if you guys prefer I can separate this question and ask another one later.~~ Lastly, here are two sample data sets upon request.               0.813010000000000   0.091178000000000   0.012728000000000     1.626000000000000   0.103120000000000   0.019204000000000     2.439000000000000   0.114060000000000   0.063494000000000     3.252000000000000   0.123130000000000   0.071107000000000     4.065000000000000   0.128540000000000   0.073293000000000     4.878000000000000   0.137040000000000   0.074329000000000     5.691100000000000   0.124660000000000   0.071992000000000     6.504099999999999   0.104480000000000   0.071463000000000     7.317100000000000   0.088040000000000   0.070336000000000     8.130099999999999   0.080532000000000   0.036453000000000     8.943100000000001   0.070902000000000   0.024649000000000     9.756100000000000   0.061444000000000   0.024397000000000     10.569000000000001   0.056583000000000   0.025222000000000     11.382000000000000   0.052836000000000   0.024576000000000     12.194999999999999   0.048727000000000   0.026598000000000     13.008000000000001   0.045870000000000   0.029321000000000     13.821000000000000   0.041454000000000   0.067300000000000     14.633999999999999   0.039596000000000   0.081800000000000     15.447000000000001   0.038365000000000   0.076443000000000     16.260000000000002   0.036425000000000   0.075912000000000      The first column is the frequencies in mHz, identical in every single data set. The second column is a good data set (good data figure one and two, panel 5, red marker) and the third column is a bad data set (bad data figure three and four, panel 5, red marker). Hope this is enough to stimulate some more enlightened discussion. Thanks you everyone.


title: How to calculate correlation between consumers according their multi-aspect and overall rating?

text: In e-commerce site, consumers would provide a overall rating with multi-aspect ratings. For example, if a consumer purchase a camera, he would provide a rating info as follows:               User    Product    Size/Wight    Appearance    Battery    Price    Display    Overall     A       Camera 1      4             2              3        5        5          4      B       Camera 2      4             2              3        5        5          2     C       Camera 3      2             4              5        1        2          4     D       Camera 4      2             4              5        1        2          2      In practice, in our dataset, a consumer would only purchase a camera. So, how can we determine the correlation between pair of users (A vs B, or A vs C, etc) using the provided heterogeneous data.


title: Identifying periods of non random behavior in time series

text: Im looking for some pointers on which topics I should be looking into in order to identify periods (of non fixed length) which deviate from randomness. I have a feeling hypothesis testing may be what I'm looking for though I haven't covered it yet. My stats knowledge is limited hence any suggested books / websites would also be appreciated.


title: Optimal basis for regression problem

text: Consider the training set $\{(x_i; y_i)\}_{i=1}^N, x_i \in \mathbb{R}^n, y_i \in \mathbb{R}$. The goal is to find regression function, like, $f(x) = \sum_{i=1}^K a_i g_i(x) + a_0$. The least-squares solution is well known, if we know $g_i$ functions. Lets consider that we have kernel regression with radial basis function (lets use gaussian), so $g_i(x) = \exp \left( - \gamma_i ||x_i - c_i ||^2 \right)$. What's is the best choice for basis function parameters $\{c_i\}_{i=1}^K$ and $ \{ \gamma_i \}_{i=1}^K $?


title: Finding patterns in data

text: I am probably looking for a definition. Imagine we have 10 variables, but we are not interested in some kind of linear relation (nor quadratic or with any curve to it). What I would like is a way to find "clusters" , patterns or combinations (whatever you want to call it). For instance, given 10 variables, let's say two (or more) people have extremely similar scores, though they are both not necessarily high nor low. I'm under the impression that this information is lost to us, while in fact this could be an interesting finding. Is there a name for trying to distinguish such interesting data patterns? Any suggestions are welcome (also for the title).


title: Why is a sample covariance matrix singular when sample size is less than number of variables?

text: Lets say I have a $p$-dimensional multivariate Gaussian distribution. And I take $n$ observations (each a $p$-vector) from this distribution and calculate the sample covariance matrix $S$. In this paper, the authors state that the sample covariance matrix calculated with $p > n$ is singular.   * How is it true or derived?   * Any explanations?


title: Understanding the two-stage choice paradigm

text: In Manski - The structure of random utility models the following example is proposed: Consider the alternative set space **_a_** = $(\alpha,\beta,\gamma)$ with the attribute representation: $X = \begin{bmatrix} & \alpha & \beta & \gamma & \cr x_1 & 1 & 1 & 2 \cr x_2 & 0.5 & -0.5 & -1 \end{bmatrix}$. The individuals/decision-makers $\sigma$ and $\tau$ are characterized by the attributes: $S = \begin{bmatrix} \sigma & \tau \cr 1 & 0 \end{bmatrix}$. Finally the function $w(x,s)$ which defines the utility of the utility maximizing agents $\sigma$ and $\tau$ is given by $w(x,s) = x_1 + x_2\cdot s$. **The aim is to compute the probabilities:** $P_t(a\in^c C)$ which is the probability that a decision-maker $t\in T$ will choose the alternative $a$ from a possible choice set $C$. By the definition of the alternative space **_a_** one can think of $C$ as any non-empty subset of any particular order of the alternatives $\alpha, \beta$ and $\gamma$. Hence there will be quite a few possible choice sets $C$ from which a decision-maker can choose. Manski denotes the choice set space, the set of all possible non-empty subsets $C$, with $\Gamma$ and $T$ represents the population of all possible decision- makers. To calculate $P_t(a\in^c C)$ additional information is given. The joint distribution $M_{\Gamma T}((C,t))$ of possible choice sets $C$ and the decision-makers $t$ is given by: $M_{\Gamma T}((\alpha,\beta,\gamma),\sigma) = \frac{2}{36}$ for each of the six ordered choice sets whose elements are $\alpha,\beta$ and $\gamma$ $M_{\Gamma T}((\alpha,\beta),\sigma) = \frac{1}{12}$ for each of the two ordered choice sets whose elements are $\alpha$ and $\beta$ $M_{\Gamma T}((\alpha,\beta,\gamma),\tau) = \frac{1}{36}$ for each of the six ordered choice sets whose elements are $\alpha,\beta$ and $\gamma$ $M_{\Gamma T}((\alpha,\beta),\tau) = \frac{2}{12}$ for each of the two ordered choice sets whose elements are $\alpha$ and $\beta$ To calculate the probabilities $P_t(a\in^c C)$ for the case that the information of $x_2$ is not observed we first need to calculate the vector $W_{Ct}$ which contains the utilities for the alternatives contained in $C$ for a decision-maker $t$. In this case this will yield: $\bar W_{(\alpha,\beta,\gamma),\sigma} = \begin{pmatrix} 1.5\\ 0.5\\ 1\end{pmatrix}$ and $\bar W_{(\alpha,\beta,\gamma),\tau} = \begin{pmatrix} 1\\ 1\\ 2\end{pmatrix}$ For the calculation of the probability $P_t(a\in^c C)$ Manski gives the following two formulas: $(1) \ P(\bar W_{Ct} | r_{Co},s_{to}) = \frac{ \sum_{ (\tilde C, \tilde t):r_{\tilde Co} = r_{Co}, s_{\tilde to} = s_{to}, W_{\tilde C\tilde t} = \bar W_{Ct} } M_{\Gamma T}((\tilde C,\tilde t))}{\sum_{ (\tilde C, \tilde t):r_{\tilde Co} = r_{Co}, s_{\tilde to} = s_{to}} M_{\Gamma T}((\tilde C,\tilde t))}$ $(2) \ P_t(a\in^c C) = \sum_{\bar W_{Ct}:\bar w_{at}\ge \bar w_{\tilde at}, \tilde a\in C}P(\bar W_{Ct} | r_{Co},s_{to})$ So clearly in order to calculate $P_t(a\in^c C)$ we need to calculate (1) first. This is done by summation over $M_{\Gamma T}((\tilde C,\tilde t))$ for all possible choice problems $(\tilde C,\tilde t)$ for particular values of $r_{\tilde Co},s_{\tilde to}$ and $W_{\tilde C\tilde t}$ namely $r_{Co},s_{to}$ and $\bar W_{Ct}$. To do that we need to consider the definition of $r_{Co}$ and $s_{to}$. Manski proposes that $X$ and $S$ are split up into observable and non-observable parts like $X = [X_o,X_u]$ and $S = [S_o,S_u]$ where the "o" donates the observable, and "u" the non observable part. Now we can write $X_o$ as $X_0 = (x_{ao}, \forall \ a\in $ **_a_** ). The definition of $r_{Co}$ is now given by $r_{Co} = (x_{ao}, \forall \ a\in C)$. So the only difference between $X_o$ and $r_{co}$ is the set on which $x_{ao}$ is defined on, since $r_{Co}$ considers all possible choice-sets $C$, $X_o$ only defines $x_{ao}$ on the alternative space **_a_**. **Now my dilemma:** Since I need to calculate (1) in order to calculate (2), I need to sum up the joint probabilities $M_{\Gamma T}((\tilde C, \tilde t))$ over possible values of $r_{Co}, s_{to}$ and $\bar W_{Ct}$. Since the value of $s_{to}$ is given by the integers $1$ or $0$ and $\bar W_{Ct}$ is just one possible utility-vector depending on which choice-set $C$ we are currently looking on, I only need to figure out what $r_{Co}$ actually is. The definition is quite clear but for this example I'm not quite sure. Since $x_2$ is not observed I thought of $r_{Co}$ as a vector of the observed outcomes of $x_1$, i.e., $r_{Co} = (1.5,0.5,1)$ for $C = (\alpha,\beta,\gamma)$. But if I consider $r_{Co}$ in this way the calculation will not work since (1) will always be 1 and the sum in (2) can get larger then 1 which is of course not possible. I'm stuck with this problem a couple of days now and any help will be much appreciated!


title: On the distribution of exam scores and the use of t-tests

text: A while back, I did some educational research where I compared exam scores between two groups. I remember reading somewhere that educational exams such as these tend to follow a normal distribution. Of course, normal distributions must necessarily be infinite and educational exam scores are not. Is it mathematically correct to say that these scores follow a "truncated normal distribution" rather than a traditional normal distribution? Does this truncation preclude one from being able to use a t-test to compare two groups?


title: Appropriate use of z test comparing two proportions

text: (Apologies for the cross-posting from math.stackexchange.com; I wasn't aware of this site until just now). I would like to test the performance of two algorithms. The second is a variant of the first. There is no prior data on the performance of either in a particular domain. The performance of both on any given sample is judged on a simple binary classification of success or failure. My null hypothesis is that there is no difference in the performance of the two algorithms. My alternative hypothesis is that the second performs better. I obtained $n$ samples and both algorithms' result on each sample was evaluated. I ended up with a proportion of success for both algorithms. Does the fact that the algorithms operated independently on the samples mean that I am now justified in using the $z$ test for comparing two proportions? Or should I have split $n$ and run a given algorithm against a portion of $n$, and the other algorithm on the other portion? Due to operational constraints the samples were/are difficult to collect, and splitting $n$ is extremely unpalatable. I also would like to be able to compare the samples directly against each other in order to get some idea of the second algorithm's limitations. What are my options? Is there some other variant of a $z$ test that applies here? Since $n$ was small to start with, should I just throw up my hands and declare this a pilot study? If I do in fact do that, am I then justified in using the $z$ test for two proportions? Any input is greatly appreciated. I'm kicking myself for not having thought this out beforehand. As you can tell, I'm not exactly an academic in this field.


title: Regression Gaussian Estimates instead of points

text: I am using a support vector regression is order to get estimates of a variable y. I want to receive a probability distribution of my estimates and not just point estimates. I want to predict Gaussian Distributions for every point where the mean should be my point and I will also give a variance. What I am thinking in order to avoid to use methods like Gaussian Processes or Probabilistic SVM is to use a history of N windows and to use the error as variance for my new estimates. For example the std of the error of the last N windows is 0.5. Then the estimate if my prediction is 5 is N(5,0.5). Is this approach right?


title: How much oversampling is justified in SMOTE algorithm?

text: I am working on an imbalanced classification problem where target to non target ratio is around 1% to 99%. Since it is an imbalanced classification problem, I have tried one strategy using SMOTE package in R. Before that additional derived variables were also added to the SAS dataset as well as some proxy business rules were added to increase the existing count of target instances. Initial data manipulations like, flooring/capping of values for few variables, variable transformations and binning of variables wrt to Gini Statistic were done prior to using the actual modeling technique Since the existing proportions of tgt to non tgt were skewed, I ran a SMOTE algorithm (perc.over = 14250% and perc.under = 100%) through R to create a balanced dataset. Post which the adjusted prior probabilities were set to the original proportions, as I have artificially inflated these proportions. Next a decision tree was on the SMOTE output to get the results. However, when I tally these results to check how many instances which were tgt = Y originally got classified, it was a disappointing number only around 15 of total tgt = Y = 1870 cases. Hence would like to take opinion on whether there is any other technique that can be used. PS: I tried using a decision tree node, embedded between start and end groups and selected boosting for 5 iterations. These results give a good True positive rate, however, the False positives are also in thousands which is again a concern. Thanks, AJ.


title: Unbiased estimator for $P(X_1=1)$

text: If $ X_1, ... ,X_n$ are IID binomial with parameters $ n$ and $p, $ find an unbiased estimator for $$G(p)=P(X_1=1)=np(1-p)^{n-1}\, .$$ I need to find this estimator so I can apply Lehmann-Scheffé theorem. I think it can't be a "crazy" estimator. But I couldn't find anything. Can I have a hint? $$\mathrm{E}[X_1(1-X_2/n)] = n p (1-p) \, .$$ Maybe it is a good way to follow.


title: Cox model with LASSO

text: Rob Tibshirani propose to use lasso with Cox regression for variable selection in his 1997 paper "The lasso method for variable selection in the Cox model" published in Statistics In Medicine 16:385. Does anyone know of any R package/function or syntax in R that does lasso with a Cox model?


title: How to compare vitamin D and glucose levels between patients and controls?

text: I want to test the hypothesis of a decreased level of vitamin D in diabetic subjects. For this I have recorded blood glucose and vitamin D levels in 40 cases and 40 controls. What kind of statistical test can I use to the above hypothesis?


title: import kernlad model to caret

text: I'd like to test the _anova rbf kernel_ included in the **kernlab** package in **caret**. Following excelent tutorial (https://topepo.github.io/caret/custom_models.html) I've come up with the following code:               SVManova <- list(type = "Regression", library = "kernlab", loop = NULL)     prmanova <- data.frame(parameter = c("C", "sigma", "degree", "eps"),                          class = rep("numeric", 4),                          label = c("Cost", "Sigma", "Degree", "Eps"))     SVManova$parameters <- prmanova         svmGridanova <- function(x, y, len = NULL) {         library(kernlab)         sigmas <- sigest(as.matrix(x), na.action = na.omit, scaled = TRUE, frac = 1)         expand.grid(sigma = mean(sigmas[-2]), epsilon = 0.000001,                     C = 2 ^(-5:len), degree = 1:2) # len = tuneLength in train         }         SVManova$grid <- svmGridanova     svmFitanova <- function(x, y, wts, param, lev, last, weights, classProbs, ...) {       ksvm(x = as.matrix(x), y = y,            kernel = "anovadot",            kpar = list(sigma = param$sigma, degree = param$degree),            C = param$C, epsilon = param$epsilon,            prob.model = classProbs,            ...) #default type = "eps-svr"     }     SVManova$fit <- svmFitanova         svmPredanova <- function(modelFit, newdata, preProc = NULL, submodels = NULL)           predict(modelFit, newdata)         SVManova$predict <- svmPredanova     svmProb <- function(modelFit, newdata, preProc = NULL, submodels = NULL)       predict(modelFit, newdata, type="probabilities")     SVManova$prob <- svmProb         svmSortanova <- function(x) x[order(x$C), ]     SVManova$sort <- svmSortanova      I then asked for the model to train some dataset:               set.seed(100) #use the same seed to train different models     svrFitanova <- train(R ~ .,                     data = trainSet,                     method = SVManova,                     preProc = c("center", "scale"),                     trControl = ctrl, tuneLength = 20,                     allowParallel = TRUE) #By default, RMSE and R2 are computed for regression (in all cases, selects the tunning and cross-val model with best value) , metric = "ROC"     #Print the results     svrFitanova      But I get the following error:               Error in train.default(x, y, weights = w, ...) :        The tuning parameter grid should have columns C, sigma, degree, eps      I don't see why this error occurs.... tune grid has four columns as requested... Any ideas? Thanks


title: How to test if some data points is drawn from a distribution with linear PDF?

text: I have some data in the range [0, 1], and from the histogram below, it seems that they might be drawn from a distribution with linear probability density function (what's the name of that kinds of distribution?). How do I estimate the parameters of that distribution and how do I test how likely the data were drawn from it in R? ![Histogram of the data](http://i.stack.imgur.com/sdUfx.png) Update: Now I know (thanks to Glen_b) if I can define the PDF (e.g., $ P(x) = ax + b$), I can test it with Kolmogorov-Smirnov test, but how do I estimate $a$ and $b$ from the sample?


title: Clinical trial question Confidence Interval

text: Data are collected in a clinical trial evaluating a new compound designed to improve wound healing in patients. The new compound is compared against a placebo. After treatment for 5 days with the new compound or placebo, the extent of wound healing is measured and the data are shown in Table. Suppose that the clinician feels that if the percent reduction in the size of the wound is greater than 50% then the treatment is a success. a. Generate a 95% Confidence Interval for the percent success in patients receiving the new compound.   b. Generate a 95% Confidence Interval for the difference in the percent success between the new compound and placbo.   c. Generate a 95% Confidence Interval for the Relative Risk of treatment success between treatments   d. Generate a 95% Confidence Interval for the Odds Ratio of treatment success between treatments.               Table:     Number of patients with percent reduction in Size of Wound     New compound n=125  None=4   1-25=11  25-50=37  51-75=32  76-100=41     Placebo      n=125  None=12  1-25=24  25-50=45  51-75=34  76-100=10      Also: a. Generate a 95% Confidence Interval for the mean age among participants assigned to the placebo.   b. Generate a Confidence interval fro the difference in mean ages in participants assigned to the experimental versus the placebo groups.   c. Generate a 95% Confidence Interval for the difference in mean BMI in participants assigned to the experimental versus the placebo groups.


title: ANOVA's and reaction time data

text: I've been reading a lot on the effect of perceptual load on selective attention but the analysis confuses me somewhat. Some of these studies flag 20-40ms differences between conditions as significant where the mean is about 500ms and my gut tells me something is fishy. Unfortunately, almost no one reports any information about the sample, anything error, SD, CI etc. Most just report sample size (usually between 10-20) and mean RT. Looking at my semi-pilot data I'd be seriously hesitant to run ANOVA's on it. It's properly positively skewed (log and inverse transformations alleviate the situation somewhat, but not always), the residuals are pretty non-normal in some cases, the dispersion of the reaction times between certain participants differs quite significantly making me wonder if I'm not going to commit a whole host of type-1 errors when running repeated -measures ANOVA's. I can see why something like fixed effects models would be preferable in these circumstance, but almost no one does that. Here is a typical example: > "A 3x2 within-subjects ANOVA was conducted on mean correct RTs and error > rates, with distractor condition (low dilution vs. high dilution vs. high > dilution/high salience) and distractor compatibility (incompatible vs. > neutral) as the within-subjects variables. All interpretations of the data > were identical using both the Sphericity assumption and the Greenhouse- > Geisser method. With respect to mean correct RTs, there was no main effect > of distractor condition (F < 1), but there was a significant main effect. of > distractor compatibility, F(1, 19) 10.06, p <.01, np2 = 0.35, indicating > that RTs were faster with neutral distractors (506 ms) than incompatible > distractors (517 ms). This main effect is qualified by a significant > distractor condition by distractor compatibility interaction, F(2, 38) = > 6.07, p <.01, np2 = 0.24." (Biggs and Gibson, 2014). Is the central limit theorem really going to make it justifiable to run ANOVA's on data from 15 participants where the distributions are positively skewed and there are differences in the variability of the reaction times between participants? This between-participant variability obviously get's chucked out the window when you aggregate these scores to produce a mean for each participant per condition. I know I'm still a noob, but I'm having a hard time making sense of this so any advice/clarification would be greatly appreciated.


title: How do you think about the central estimate when the confidence interval is asymmetric?

text: I'm using a wild bootstrap to create confidence intervals around fitted values of the following model, for a specific combination of the factors, as x varies across its range.               mod = gam(log(lhs)~ s(re,bs="re",by=dum)         +   as.factor(fac1)         +   as.factor(fac2)         +   as.factor(fac3)         +   s(x,by=as.factor(fac1))         +   s(x,by=as.factor(fac2))         ,method="REML"         ,data = dat)      I am not simply doing `predict.gam(mod, newdata=whatever,se.fit=TRUE)` because the dependent variable is a log, and the SE*1.96 overlaps zero. What I want is a distribution of $lhs^{-1}|x$ that I can use in a monte carlo simulation, where I'd specify x and take a draw from the conditional distributon of $lhs$. So I am using a wild bootstrap -- taking the residual vector, multiplying it by the Mammen distribution, adding it to the fitted values, and refitting the model to the resultant $lhs^*$. To be precise, I am using a cluster wild bootstrap, because my data is clustered. When I do so, I get this: ![enter image description here](http://i.stack.imgur.com/dBxlZ.png) This is 50-odd bootstraped fitted values. I want to sample from this when doing the monte carlo -- giving an $x$ and getting a possible value of $lhs|x$. **My question: since the central estimate is clearly not at the center of the distribution of fitted values, what, if any, information or value does the central estimate contain? Or, should I take this asymmetry as an indication that I am either using the wrong tool, or using it improperly?** If the latter, I'd greatly appreciate an explanation of why what I'm doing is wrong. (Note: this is a plot of $lhs^{-1}|x$. I have also looked at plots of $log(lhs)|x$, which is what the model specifies. It is are similarly non- symmetric, so my problem doesn't appear to be that the CI is simply symmetrical under transformation.) (One more note: the distribution shown -- above the central estimate on the left and below it on the right - is highly consistent with my priors). Advice very much appreciated.


title: Prediction with Bayesian networks in R

text: I've been trying to teach myself about Network Analysis, and I've been able to develop DAG charts in R. However, I've looked through three or four R packages and have seen little in the way to a function to generate joint probabilities for the network. The DAG plot tells me about the variables in relation to one another, but I'm more curious about the probabilities and haven't found a way to do that in R. If anything, there seem to be a lot of packages dedicated to generating plots or focused on inference, and I want to know how I can get the probabilities for the networks.               library("bnlearn")     library("Rgraphviz")          dat=data.frame(won=c(1,0,0,1,0,0), sold=c(0,0,0,1,0,0), insured=c(0,0,1,0,0,1),                     credit=c("POOR","FAIR","GOOD","FAIR","FAIR","GOOD"))     dat$won = factor(dat$won)     dat$sold = factor(dat$sold)     dat$insured = factor(dat$insured)     dat$credit = factor(dat$credit)           highlight.opts <- list(nodes = c("won","sold","insured","credit"),                            col = "red", fill = "grey")     bn.hc <- hc(dat, score = "aic")     graphviz.plot(bn.hc, highlight=highlight.opts)


title: Is it possible to specify a custom likelihood function using the "link" variable in the `glm` package?

text: I am trying to create a custom likelihood function that uses different likelihood functions depending on whether an observation is endogenous or exogenous. Is it possible to hardcode my function using the link variable in the `glm` function in the `glm` package? Or would I have to do everything by hand and then use an optimization procedure like `optim`? Thank you!


title: Post hoc test of adjusted means (after ANCOVA)

text: When doing ANCOVAs, software usually gives both unadjusted means as well as adjusted means. However, the post hoc tests for most statistical programs are limited to unadjusted means. Yet, sometimes, if one is doing an ANCOVA, post hoc comparisons of adjusted means is what is desired, not the unadjusted comparisons. Or am I mistaken – should post hoc tests of adjusted means be avoided for some reason? Post hoc tests for adjusted means have been around for quite some time, such as the GT2 method of Hochberg (1974). A modified version of the GT2 test appears in Sokal and Rohlf's Biometry (1981 2nd edition, and presumably in newer versions) specifically for adjusted means. There is also the T-method (Spjotvoll and Stoline 1973), also modified in Sokal and Rohlf (1981) for comparing adjusted means. Why then do statistical programs not give post hoc tests for adjusted means, if one wants to know which treatment (A vs B, A vs C, etc) differ, after removing some of the variance of the covariate? Hochberg, Y. 1974. Some Conservative Generalizations of the T-Method in Simultaneous Inference. Journal of Multivariate Analysis, 4:224 -234. Spjotvoll, E, Stoline, MR. 1973. An extension of the T-method of multiple comparison to include the cases with unequal sample sizes. J Am Stat Assoc. 68:976-978.


title: Cross validation as a measure of SVM generalizability

text: I am trying to implement the following paper: Learning the Kernel Matrix with Semidefinite Programming using cvx toolbox for matlab. My question is about how exactly I can use this method with cross-validation to generate some measure of average accuracy or "generalizability". My situation is the following: I have 100 trials of some task from a person. These are all labeled, in that we know what type of trial we gave the person to perform: `A` or `B`. I now want to build a classifier to see if we can tell them apart without the labels, in the hopes of generalizing into unlabeled data in the future. So, I went about implementing the 2-norm soft margin equation in the paper (equation 41) using cvx to learn the kernel matrix to try and learn the ideal kernel matrix. I am unclear how exactly this fits in with cross-validation. Does it make sense to split my data into train/test, then learn the kernel matrix using the technique described and test on my test data, and then shuffle my data again into a new train/test set and re-learn the kernel? How does this speak towards the generalization of this when Im essentially learning a new kernel every time? Imagine I do a 5 fold cross validation and get something like: 50, 75, 85, 95, 100%? What can I conclude from that? That my classifier has 81% accuracy? What about the 50% and 75% trials? What if my trails have some mean, but a HUGE variance? How can I be confident? Also, which one of the five kernel matrices that I've now learned, do I use on the newly collected unseen data?


title: Exploratory data analysis for discrete data

text: I am using a probit and a logit model for obtaining the choice probabilities of some data. What kind of plots can be useful to conduct a exploratory data analysis for these data? Here is the simulated data:               set.seed(1)     x.first <- rbinom(1000,1,0.5)     x.second <- rbinom(1000,1,0.5)     y <- rbinom(1000,1,0.5)     delta.x <- x.first-x.second      and my models are                logit(p) <- delta.x beta           phi^{-1}(p) <- delta.x beta      where phi^{-1} is the inverse of the standard normal distribution `cdf`and beta is the regression coefficients.


title: Why should we care about rapid mixing in MCMC chains?

text: When working with markov chain monte carlo to draw inference, we need a chain that mixes rapidly, i.e. moves throughly the support of the posterior distribution rapidly. But I don't understand why we need this property, because from what I understand, the accepted candiate draws should and will concentrated in the high density part of the posterior distribution. If what I understand is true, then do we still want the chain to move through the support ( which includes the low density part) ? In addition, if I am using MCMC to do optimization, do I still need to care about rapid mixing and why? Thank you for share your thoughts!


title: Multiple linear regression with restricted ranges in R

text: I'm trying to do a simple least-squares fit on some data. The formula is simply: $Y = C_0X_1 + C_0X_2 + C_0X_3 + C_0X_4 + C_0X_5$ I have 24 rows of $Y$'s and $X$'s, and I'm trying to fit the $C$'s. The problem, however, is that I need to constrain the ranges of the $C$'s to be non-zero, regardless of the quality of the fit. I've read the docs on `lm` in R, but I cannot find the option to constrain the variable ranges. Any ideas?


title: Good Tutorial for Restricted Boltzmann Machine (RBM)

text: I’m studying the Restricted Boltzmann Machine (RBM) and having some issues understanding log likelihood calculations w.r.t. parameters of the RBM. Even though a lot of research papers on RBM have been published, there are no detailed steps of the derivatives. After searching online I was able to find the following document. http://image.diku.dk/igel/paper/AItRBM-proof.pdf However, the details of this document are too advanced for me. Can somebody point out me to a good tutorial/lecture note about RBM? Edit: @ David Confusing sections are shown below ( equation 29 in page 26). ![log likelihood calculations](http://i.stack.imgur.com/xaIcq.png)


title: understanding variation of data and visualization

text: 1) What is the threshold or measure to say that the variation of data is high? Typically we say high or low variation by comparing with another data set. But there is only one data set. 2) We can use boxplot to visualize standard deviation, mean, etc., of the data. However, I am curious to know whether there is any visualization method or distribution method in R to _stress_ that the variation is high in a graph. For example, I have a data set of 64 points and it is one dimensional. Please see data below. Assume the values are always positive. Its mean is 32.8, standard deviation(sd) is 37.7, and coefficient of variation(cv) is 1.2. May be these are silly questions. But I won't get peaceful mind if I don't get a clarification. Please help me. Thank you.               0 19  0 83 88 97  0 50  0 82 89  0  0  0  0  0  0 17 33 34  0  0 50 59 84 0 67 83 89  0 17  0 66  1 59 89  0 14 68  0  0  0 95  0  0 89  4 85 89 25 50 83 34  0  0  0  0  0  0  0 96  0 89 19


title: The expressive power of PoE over MoE

text: I am reading Hinton's paper "Training Products of Experts by Minimizing Contrastive Divergence" about training PoE. There is an explanation of PoE's expressiveness, which stats "because the posterior distribution cannot be sharper than the individual models in the mixture and the individual models must be broadly tuned to allow them to cover the 35-dimensional manifold". I can understand the first half but I don't quite understand the second half: "the individual models must be broadly tuned to allow them to cover the 35-dimensional manifold". Is it possible that one individual model in MoE is sharp at some point while other models be small at the same point? If it is possible, this will make the posterior be very sharp at that point.


title: How to combine several categorical variables into one quality of health measure?

text: I am new to SPSS, and am trying to use SPSS to generate a variable on the quality of health service available to the residents of an area. Below are the categorical variables that could tell me the quality of health available to them.   * no insurance covering health care cost    * no access to ambulance service    * poor attention by medical personnel    * poor medical equipment    * poor drug provision    * time taken from dwelling to health centre more than 30mins    * poor waiting time at health facility    * health facility not available within neighbourhood  How can I generate a number value on the quality of health of this population from the number values of these variables, bearing in mind that several respondents may affirm each variable as true.


title: Are the $1-SSe/SSt$ and $cor^2$ calculations of $R^2$ always equivalent?

text: I am trying to calculate the $R^2$ value for a production constrained spatial interaction model, using Fotheringham and O'Kelly (1989) as my guide. I get dramatically different values for R-Square, depending on whether I calculate it as `r-square <\- 1 - SSe/SSt` or `r-square <\- cor(x, y)^2`. Is this result expected? Of course, I may well be miscalculating this somewhere along the line. I want to use r-square as a (flawed but nevertheless useful and widely understood) measure of goodness of fit, as recommended by Fotheringham & Knudsen (1987). A reproducible example is below. I've saved my model output to a csv, to save space here.               predobs <- read.csv("http://dl.dropbox.com/u/66606821/pred_obs.csv")     sst <- sum((predobs$obs - mean(predobs$obs))^2)     sse <- sum((predobs$obs - predobs$pred)^2)     (r.square.1 <- 1 - (sse/sst))     (r.square.2 <- cor(predobs$obs, predobs$pred)^2)


title: What are chunk tests?

text: In answer to a question on model selection in the presence of multicollinearity, Frank Harrell suggested: > Put all variables in the model but do not test for the effect of one > variable adjusted for the effects of competing variables... Chunk tests of > competing variables are powerful because collinear variables join forces in > the overall multiple degree of freedom association test, instead of > competing against each other as when you test variables individually. What are _chunk tests_? Can you give an example of their application in `r`?


title: Improving an unbiased estimator using the Rao-Blackwell theorem

text: Given data following a Bernoulli distribution $Y_1,...Y_n \sim B(1,p)$, we want to measure $\theta = \text{Var}(Y_1)$ I found this unbiased estimator $S^2 = \frac{\sum{(y_i-\bar{y}})^2}{n-1}$ I need to improve this estimator using the Rao-Blackwell theorem and determine if the new estimator is UMVUE for $\theta$. I started with identifying the unbiased estimator sufficient statistics. $$T(y)= \frac{\sum{(y_i-\bar{y}})^2}{n-1}$$ $$w(y)= \sum{(y_i-\bar{y}})^2$$ Therefore $$T_1=E(T/W)=E(y_1/\sum{(y_i-\bar{y}})^2)$$ I'm not sure regarding the sufficient statistic ($w(y)$). Is it okay or am I supposed to just use $\sum{y_i}$ Next I tried to apply the Rao-Blackwell theorem, but I'm not sure about it at all $$y_1/T(y)\sim B(1,P(y_1=1,\sum{(y_i-\bar{y})}^2=t))$$ $$\frac{P(y_1=1,\sum{(y_i-\bar{y})}^2=t)}{P( \sum{(y_i-\bar{y}})^2=t)} = \frac{P(y_1=1)\cdot P(\sum{(y_i-\bar{y}})^2=t-1)}{\sigma^2 \cdot \chi^2_{n-1} }$$ I apologize if this is rubbish, it's just that I'm trying to figure out exactly how to go about this from what we were taught (or not, for that matter)


title: Advice on how to analyse “customer-data” in R

text: consider the following example data:               df1 <- data.frame(customer=c(rep("customer1",5),rep("customer2",10),rep("customer3",7)),                       money_spent=sample(22))          df2 <- data.frame(customer=c("customer1","customer2","customer3"),                       origin=c("US","US","UK"),                       industry_sector=c("IS1","IS2","IS3"),                       currency=c("USD","USD","GBP"))      My actual data consists of about 200000 rows and I would like to examine it in terms of, for instance, do customers from the US spent more money compared to customers from other countries. I would also like to see whether the amount of money spent depends on the industry sector and so on. I have some more explanatory variables apart from origin, industry sector and currency which I would like to look into. Also the number of records for the customers differ so that it might make sense to average the money spent for each customer. I am not sure about how to best analyse this data. I first thought of cluster analysis, in particular, hierarchical clustering but am not sure whether it can be applied to such data and, in particular, how to structure the data to be put into the function. The R function `hclust` takes a matrix as the input but how would I structure such a matrix in terms of my data? Could k-means clustering be a better alternative? Another approach would probably to analyse this data using boxplots and an one-way ANOVA approach to see whether "money spent" differs between different countries or industry sectors. However, this approach does not test whether the variables are dependend on each other. To look into this, I have been advised to apply a decision tree first and then do some statistical significance analysis. However, from what I have read so far I cannot see how decision trees can help me to detect variable dependencies. So, I am wondering whether there any other/better techniques/functions out there which are more suitable for such data? Maybe a time series analysis is more appropriate since we have also recorded the dates when customers spent money.


title: Minimum points required for fitting curve

text: If I need to build yield curve, n what is the minimum number of data points necessary for any interpolation method, e.g. say cubic spline or Nelson–Siegel, etc. I ma not sure if I can raise this query here, but if someone can guide, it will be a great help. Thanks.


title: Getting the variance-covariance matrix of regression coefficients in GEE

text: I fitted a GEE model using the function `genZcor` with user defined correlation matrix. I want to get the var-cov matrix of the regression coefficients. But the output provides only limited information. I would be very much thankful if you could kindly let me know how to get it since I am struggling lot getting this.


title: Degrees of freedom of a linear function of four independent random variables

text: This is basically a modified version of a previously unanswered question I posed here, as I'm starting to have second thoughts about my methodology. In summary, I have four random variables: Xideal, Xnonideal, Yideal, Ynonideal, which I have sampled a different number of times (> 30). I am able to get the difference between Xideal and Xnonideal (denoted `X`) as well as the difference between Yideal and Ynonideal (denoted `Y`), which I can display as mean difference and standard error:                    N_ideal  N_nonideal  Mean Difference  Standard Error     X    64       62          23.5             1.56     Y    71       67          27.1             1.67      I'm now looking to get the p-value to assess the difference between X and Y. The t-statistic if I'm not mistaken is simply $$ t_s = (27.1 - 23.5)/\sqrt{1.56^2 + 1.67^2}; $$ I guess the question is now, how do I calculate the degrees of freedom (df)? Edit: I calculated the standard error with the following formula $$ ste_{X} = \sqrt{ \sigma^2_{X_{ideal}}/N_{X_{ideal}}+\sigma^2_{X_{nonideal}}/N_{X_{nonideal}}} $$ Edit 2: Ok, to sum this up, I have four independent random variables, I need the df of $$ Z = (X_{i}-X_{n}) - (Y_{i}-Y_{n}) $$ Looking at the Welch–Satterthwaite equation, would my answer be: (if $$Nt_x = N_x - 1$$) $$ v_Z = \frac{(Var(X_{i})+Var(X_{n}) + Var(Y_{i})+ Var(Y_{n}))^2}{(Var(X_{i}))^2/Nt_{Xi}+(Var(X_{n}))^2/Nt_{Xn}+(Var(Y_{i}))^2/Nt_{Yi}+(Var(Y_{n}))^2/Nt_{Yn}} = 236 $$ or $$ v_Z = \frac{(Var(X_{i})-Var(X_{n}) - Var(Y_{i} )+ Var(Y_{n}))^2}{(Var(X_{i}))^2/Nt_{Xi}+(Var(X_{n}))^2/Nt_{Xn}+(Var(Y_{i}))^2/Nt_{Yi}+(Var(Y_{n}))^2/Nt_{Yn}} = 16 $$


title: which algorithm uses a multidimension array of average probabilities

text: **Is there an algorithm that returns the mean probability of response from a multi-dimensional matrix?** For example, if I have a set of features:   * customerClass = `poor` / `middle` / `rich`   * category = `11` / `20` / `35`   * offerValue = `$1` / `$2` / `$3` and a response:   * usesCoupon = FALSE / TRUE Assuming I have enough samples, I want to predict the probability for a new sample (customerClass = `rich`, category = `20`, offerValue = `$1`) as the average probability for the population with the same three-dimensional value `(rich, 20, $1)`.


title: Generating gamma random field with given covariance matrix

text: I have to generate multivariate gamma distributions with given positive- definite covariance matrix. Anyone can suggest me a method? Thanks,


title: Interpretation of the range parameter in a Gaussian Process

text: I have been teaching myself about Gaussian Process (GP) modeling for a while now and although it is "easy" to estimate the range parameter (sometimes called the length-scale) in the GP I am actually trying to gain a physical understanding/intuition for what the range parameter actually represents. Just so that my question is clear, if our GP had the following correlation function: $$R(\phi) = \exp\left(-\frac{|x_i-x_j|}{\phi}\right)$$ then $\phi$ would be the range parameter. I am also interested in gaining intuition on knowing whether or not the value of $\phi$ that I estimate is reasonable or not given my data and true function that I am trying to model.


title: Growth models for a population

text: I'm studying the growth of a population (users of a website). I have the user count for each time block (which is 2 weeks long). Now, I'd like to understand if this growth follows a well-known growth curve. To start with, I ran this regression on R on 30 observations with `glm(tot_users ~ block_id,family=binomial(logit), data = df)`. These are the observations: > > 25,83,111,164,251,370,557,815,1154,1513,2032,2605,3590,4904,5718,6602,7628,8727,9471,10263,11047,11799,12441,13040,13634,14168,14582,15143,15649,16164,16472 This is what I get:               Coefficients:     (Intercept)   block_id10   block_id11   block_id12   block_id13   block_id14   block_id15   block_id16   block_id17   block_id18   block_id19    block_id2   block_id20   block_id21   block_id22   block_id23   block_id24   block_id25   block_id26   block_id27   block_id28   block_id29    block_id3   block_id30   block_id31    block_id4    block_id5    block_id6    block_id7    block_id8        -9.397e+00    6.331e+00    6.667e+00    7.148e+00    7.398e+00    7.714e+00    7.968e+00    8.151e+00    8.364e+00    8.623e+00    8.979e+00    1.958e-13    9.263e+00    9.430e+00    9.591e+00    9.817e+00    9.985e+00    1.020e+01    1.048e+01    1.067e+01    1.098e+01    1.126e+01    6.932e-01    1.196e+01    3.196e+01    2.946e+00    3.468e+00    4.754e+00    5.555e+00    5.896e+00         block_id9         6.150e+00            Degrees of Freedom: 30 Total (i.e. Null);  0 Residual     Null Deviance:      17.57      Residual Deviance: 3.167e-10    AIC: 74.87      If I plot this model, I get a surprisingly close approximation. The residual deviance seems incredibly low, indicating that the model is very good (is it?). My questions are: 1) Can I assume that my variable growth follows a logistic curve? 2) What other curves do you think I should try to fit? 3) What is the usual validation process to make the regression results publishable? What measures are usually reported? Thanks for any hint. Mulone


title: Results of bootstrap reliable?

text: I am using the bootstrap algorithm to compute standard errors of the estimates of my normalmixEM output. I am not really sure if they are reliable? My code is (data here):               # load package     install.packages("mixtools")     library(mixtools)               B = 1000 # Number of bootstrap samples     mu1sample <- mu2sample <- sigma1sample <- sigma2sample <- lambdasample <- vector()          # Bootstrap     for(i in 1:B){       print(i)       subsample = sample(mydatatest,rep=T)       normalmix <- normalmixEM(subsample, mu=c(-0.002294,0.002866),sigma=c(0.00836,0.02196), lambda=c(0.6746903,(1-0.6746903)),k=2, fast=FALSE, maxit=10000, epsilon = 1e-16, maxrestarts=1000)       mu1sample[i]    = normalmix$mu[1]      # $       mu2sample[i]    = normalmix$mu[2]      # $       sigma1sample[i] = normalmix$sigma[1]   # $       sigma2sample[i] = normalmix$sigma[2]   # $       lambdasample[i] = normalmix$lambda[1]  # $     }     # standard errors          sd(mu1sample)     sd(mu2sample)     sd(sigma1sample)     sd(sigma2sample)     sd(lambdasample)          # show distribution of the bootstrap samples     hist(mu1sample)     hist(mu2sample)     hist(sigma1sample)     hist(sigma2sample)     hist(lambdasample)      This gives the following pictures: mu1 ![m1](http://i.stack.imgur.com/L0n8F.png) mu2 ![m2](http://i.stack.imgur.com/RG2OA.png) sigma1 ![s1](http://i.stack.imgur.com/UnUcm.png) sigma2 ![s2](http://i.stack.imgur.com/iipWM.png) lambda ![la](http://i.stack.imgur.com/kHYnM.png) EDIT: If you look at my variable, the mydatatest and use a KD to show the distribution with the following code               plot(density(mydatatest),col="red",main="",lwd=2,cex.axis=1.2,cex.lab=1.2)      it looks like ![tt](http://i.stack.imgur.com/kaAYE.png) 2nd EDIT: I now included the mus and sigmas to be fixed. I updated the code and the pictures. Now again my question, what do you think about it?


title: Variance and covariance of binary data

text: The variance of a set of $n$ binary variables $D = <x_1, \ldots, x_n>$ is $$ {\rm var(D)} = \frac{k(n-k)}{n^2}, $$ where $k$ denotes the number of $1$s in $D$ (see http://capone.mtsu.edu/dwalsh/VBOUND2.pdf). With Matlab the variance of $<1, 0, 0, 1, 0, 0>$ is 0.2667, while using the above formula the variance is 0.22. What is the reason for that? Moreover, what is the simplified version of covariance formula between two binary variables?


title: Implementation of Krippendorff's α for unitizing continuous data

text: Krippendorff has a paper (free and official versions) describing an adaptation of his measure of reliability α to continuous data, such as text. Specifically, it calculates the reliability of both unitizing (i.e., how the text is split into coded and uncoded units) and coding (i.e., which codes are assigned to what units). I'm wondering if this has been implemented, say, in R. In his paper, Krippendorff mentions that a "computer program for calculating these [adaptations of α] is currently being developed." However, I could find no mention of or link to such software, neither on Krippendorff's site nor elsewhere. The paper is detailed enough that it would allow me to implement it myself, but if it's already been done no sense in reinventing the wheel. Thanks in advance.


title: Why should all Cross-Validation results be higher than the result on the test dataset?

text: Sorry, I'm not an expert and my question could be fundamentally wrong. I've read this interesting question because I also was wondering whether to train the model again after cross-validation. Now, after boosting regression trees with this library over 4/5 of the training set, I'm doing a 5-fold CV against a particular objective function. After re-training on the whole training dataset, I use that final model to predict the binary classification of the test dataset. I'm a bit surprised to see that the evaluation metric for the test set is significantly worse than expected from the validation results. Do you have any suggestions about what could be going wrong? My second (maybe related) question is why does increasing the rounds of boosting result in worse predictions?


title: A Kalman Filter for estimating z-scores?

text: I have been struggling to fit the following problem into a linear state space model for a Kalman Filter (KF). I'm having a hard time seeing what I'm doing wrong. I suspect I'm violating some law of KF state space modeling. Perhaps something is non-gaussian, non-linear, non-kosher? I have data from players who score points in a game each day. The game involves both skill and luck. With a KF in mind I am thinking about skill as the unobserved state and luck as part of the observation noise. Here is some (contrived) data:                       Points Scored     Player  Day 1   Day 2   Day 3     1       10      118     50     2       15      100     49     3       20      130     51      I am also assuming players are endowed with a skill level (relative to the other players) which does not change over time. Further, points scored are exogenously and randomly mean shifted/scaled each day. As you can see in the example data above, Player 3 is more skilled because he scores more points than the other players. Also, the distribution of points moves and scales -- that is to say means and variances of points scored change from day to day. I'll index players $i=1,\dots,N$ and time $t=1,\dots,T$. Let $y_{it}$ be the points scored by player $i$ at time $t$. My first approach was to estimate skill using something like a z-score (or a t-stat). For each player at each time step I compute: $z_{it} = \frac{y_{it} - \bar{y_t}}{s_{y_t}}$ Then I average across time to estimate skill: $\widehat{skill_i} = \bar{z_i}$ As far as I can tell, so far so good. Now I want to Kalmanize this thing so I can update estimates of skill as new data arrive. I start to think about the state and its relationship with the observable points scored. But this is where I get stuck. Suppose the true, unobservable skill level for a player $i$ is $\theta_i$ and this comes from a gaussian distribution with mean and variance $(\mu_\theta, \sigma_\theta^2)$. Luck is $\epsilon_{it}$ with mean and variance $(0,\sigma_{\epsilon}^2)$ Further suppose the common daily mean shift and scaling factor are $\alpha_t$ and $\gamma_t$, respectively, where both come from gaussian distributions with means and variances $(\mu_{\alpha}, \sigma_{\alpha}^2)$, and $(\mu_{\gamma}, \sigma_{\gamma}^2)$. I model the observed points scored by player $i$ at time $t$ as: $y_{it} = \alpha_t + \gamma_t (\theta_i + \epsilon_{it})$ The state transition equations are easy: $\theta_{it} = \theta_i$. Each player's skill is assumed to be constant across time. The player's luck is the observation noise. Although it seems like I have everything defined right, I can't quite seem to figure out how to make $\alpha_t$ and $\gamma_t$ fit into a KF model - should they be part of the state? If so, how can I get a linear KF to multiply states? If they aren't part of the state, but are time-varying coefficients, how can I make that work? And of course I should also be asking: Is the Kalman Filter even an appropriate estimation strategy for skill? Help! So many thanks in advance!


title: Post hoc $\chi^2$ test with R

text: This is a quite simple question but I don't find any good, clear, precise answers: I'm looking for a way to perform _post hoc_ test on a chi$^2$ test. I have 2 variables : var1 : good/fair/poor and var2: a/b/c This is the contingency table :                     a   b  c     good 120  70 13     fair 230 130 26     poor  84  83 18      with R :               mat <- matrix(c(120,230,84,70,130,83,13,26,18),3)     dimnames(mat) <- list(c("good","fair","poor"),c("a","b","c"))     mat      I perform $\chi^2$               chisq.test(mat)      I obtain p =0.022, so there is a relation between var1 and var2... but which one? How can I handle this? (contribution analysis?, pairwise prop test?)


title: Overfitting on Missing Value Imputations

text: When performing Missing Value imputations, should we be concerned about overfitting the data? Why or why not? For example: If I impute a variables missing value using a CART regression tree, should I be concerned about my tree being too complex and overfitting that particular data? When it comes time to use testing data, should I use that same tree model? Or should I build a new tree based off of my testing dataset?


title: How to measure the similarity of k-means clustering using different datasets?

text: I run k-means clustering on my dataset (100 samples in total) and partition the data into k=5 clusters. Then I want to test how robust of the k-means can be; however, I haven't got more new data samples. My idea is:   * Take the first sample out and run k-means on the rest of 99 samples.   * Loop over the step described above for each sample (e.g., take out the 2nd sample at the 2nd iteration), and run the k-means 100 times in total. My question is how to measure the similarity of the 100 k-means results? I am thinking of get the statistics of silhouette coefficients. Does that make sense? Thanks. A.


title: Conditional or unconditional exact test in R

text: I have a 2x2 contingency table and i want to calculate if the pair inside is significantly different. i made a matrix like the following named raw_matrix                         CNS random     Not_H3K4  343  28825     H3K4      11   2014      Create this matrix , thus:               raw_matrix = structure(c(343, 11, 28825, 2014),          .Dim = c(2L, 2L), .Dimnames = list(         c("NotH3K", "H3K"), c("CNS", "Random")))      as i searched, unconditional exact test like Barnard’s and Boschloo’s exact tests are the most powerful test for this end. i installed the 'Exact' package and tried to do the test using this command:               exact.test(raw_matrix)      it took more than half an hour on a 64GB ram and 3.5 GH CPU computer and finally it gave the following error:                   Error: cannot allocate vector of size 42.0 Gb     In addition: Warning messages:     1: In matrix(A[xTbls + 1, ] * B[yTbls + 1, ], ncol = length(int)) :       Reached total allocation of 61417Mb: see help(memory.size)     2: In matrix(A[xTbls + 1, ] * B[yTbls + 1, ], ncol = length(int)) :       Reached total allocation of 61417Mb: see help(memory.size)     3: In matrix(A[xTbls + 1, ] * B[yTbls + 1, ], ncol = length(int)) :       Reached total allocation of 61417Mb: see help(memory.size)     4: In matrix(A[xTbls + 1, ] * B[yTbls + 1, ], ncol = length(int)) :       Reached total allocation of 61417Mb: see help(memory.size)      then i installed 'Exact2x2' package and did the test using this command:               exact2x2(raw_matrix)      which gave me the following results:                   Two-sided Fisher's Exact Test (usual method using minimum likelihood)          data:  raw_matrix     p-value = 0.006433     alternative hypothesis: true odds ratio is not equal to 1     95 percent confidence interval:      1.2028 4.2424     sample estimates:     odds ratio        2.178631       but as i read in the 'Exact'package tutorial , the fisher exact test which is a conditional exact test is not so powerful. finally i did the normal chi square test using the command chisq.test(raw.matrix) which gave the following results that is different from fisher test's results:                   Pearson's Chi-squared test with Yates' continuity correction          data:  test_1     X-squared = 6.2045, df = 1, p-value = 0.01274      im a Geneticist and not an expert in statistics, i appreciate if anybody could tell me what is the best strategy here to do this test


title: Choice of dependent variable: Differencing or Controlling?

text: I was running some analysis where I suspect that a treatment $D$, has opposite effects on two variables $Y^A$ and $Y^B$. To show that, I was thinking about two strategies: **1. Differencing** Running the regression using the difference between $Y^A$ and $Y^B$ as my dependent variable, namely: $(Y^A_i-Y^B_i)=\delta^d D_i+\epsilon_i, $ where $\delta^d$ is the coefficient of interest. **2. Controlling** Alternatively, I might regress on $Y^A$ while controlling $Y^B$ and the interaction between $D$ and $Y^B$ (and possibly higher order terms of $Y^B$ for nonlinearity), which is something like: $Y^A_i=\delta^cD_i+\beta_1Y^B_i+\beta_2(D_i\times Y^B_i)+\cdots(\textrm{higher order interactions})\cdots+\epsilon_i$. In my data, $Y^A$ and $Y^B$ are correlated ($\rho\simeq 0.2$).In both cases, the quantity that I try to get at is how much extra change in $Y^A$ does the treatment in addition to the treatment's influence on $Y^B$ (which in turn influences A). Could someone enlighten me about the statistical properties of each approach, as well as which one is appropriate to use for my purpose here? Or, should I be using a different (third) method? Thanks in advance!


title: How to apply predict functionality of R in SAS?

text: do we have any proc or way that we can apply a model built on train dataset to test dataset using the predict function similar to R in SAS?


title: Can PCA scores be used as dependent variable?

text: I am working on a research project where I have several questions from a survey data that measures the same underlying quantity (my dv), possibly each with some measurement error. I was thinking about using PCA method to extract the most relevant component. My initial analysis turns out to give what I wanted, but I'm not sure if this is a generally accepted practice in statistics. Could someone let me know if there is any problem with this, and if there is any publication that has taken a similar approach? Thanks a lot


title: Is the Poisson distribution stable and are there inversion formulas for the MGF?

text: First, I have a question about whether the Poisson distribution is "stable" or not. Very naively (and I'm not too sure about "stable" distributions), I worked out the distribution of a linear combination of Poisson distributed R.V.'s, using the product of the MGF. It looks like I get another Poisson, with parameter equal to the linear combination of the parameters of the individual R.V.'s. So I conclude that Poisson is "stable". What am I missing? Second, are there inversion formulas for the MGF like there are for the characteristic function?


title: Inequality convolution distribution function

text: $F$ is the distribution function of a positive random variable. Why does the following hold for every $x \in (0,\infty)$: $F^{*2}(x) \leq F^2(x)$, where $F^{*2}=F*F=\int\limits_0^x F(x-y)dF(y)$


title: How important is it to define the nested design in a within subject ANOVA?

text: I am asking to find out why a factor may become non significant after I define the nested design. The following are the two aov analysis in R. Obs (subject #) is the within subject factor and stim3 & label (1st & 2n stimuli manipulations) are the nested within subject. stim3 & label are nested within obs. The analysis shows that 'label' is only significant when the analysis is done without defining the nested design. I wanted to ask why this is the case and how sound is the choice to define the nested design.                > summary(aov(resp2 ~ stim3*label+Error(obs/(stim3*label)), tl ) )              Error: obs                   Df Sum Sq Mean Sq F value Pr(>F)         Residuals 11   5.25   0.478                             Error: obs:stim3                   Df Sum Sq Mean Sq F value  Pr(>F)             stim3      1  165.4   165.4     425 3.8e-10 ***         Residuals 11    4.3     0.4                             ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               Error: obs:label                   Df Sum Sq Mean Sq F value Pr(>F)           label      1   6.77    6.77    4.77  0.051 .         Residuals 11  15.61    1.42                          ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               Error: obs:stim3:label                     Df Sum Sq Mean Sq F value Pr(>F)         stim3:label  1  0.026  0.0255    0.15   0.71         Residuals   11  1.893  0.1721                             Error: Within                     Df Sum Sq Mean Sq F value Pr(>F)         Residuals 1872    279   0.149                             > summary(aov(resp2 ~ stim3*label+Error(obs), tl ) )              Error: obs                   Df Sum Sq Mean Sq F value Pr(>F)         Residuals 11   5.25   0.478                             Error: Within                       Df Sum Sq Mean Sq F value  Pr(>F)             stim3          1  165.4   165.4 1048.60 < 2e-16 ***         label          1    6.8     6.8   42.95 7.2e-11 ***         stim3:label    1    0.0     0.0    0.16    0.69             Residuals   1905  300.4     0.2                             ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       Thanks for reading!


title: Mixed model in SPSS with random effect and repeated measures

text: I am working on analyzing a dataset that involves repeated measures data. The data was previously analyzed by a colleague using custom code written in C++, but I have expanded the dataset and am trying to analyze it in SPSS in a similar way. The dataset has a subjects variable that I want to specify as a random-effects variable and two within-subjects variables with two levels each. From what I've read so far, the Mixed Model command in SPSS seems to be the most appropriate way to analyze this data. I'm unfamiliar with analyzing the data in this way but I've tried to read up on it and I will post the syntax I used to supplement my questions. The experimental design is like this: I have one variable, subjects, that I want to include as a random factor. The subjects saw images followed by letter strings on a computer screen and had to decide whether the string was a word or non-word. I have two within-subjects factors: "wordness" (whether the string was a word or non-word) and "relatedness" (whether the string was related to the image or unrelated). Subjects will perform multiple trials, so they will have responses for all four combinations of wordness and relatedness: related word, related non-word, unrelated word, and unrelated non-word. Together, I think, this forms a nested repeated measures design. The dependent variable is reaction time. My biggest issue is with how to specify "subjects" as a random effects variable. I've tried adding it as a subjects variable through the GUI (yielding the following syntax): MIXED Rt2Adj BY Relatedness Wordness Subjects   /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1)   SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)   /FIXED=Relatedness Wordness Relatedness*Wordness | SSTYPE(3)   /METHOD=REML   /PRINT=SOLUTION   /REPEATED=Relatedness*Wordness | SUBJECT(Subjects) COVTYPE(UN). This runs without errors, but I notice that there are no random effects variables listed in the Model Dimension output table. If I add subjects as a random effects variable in addition to specifying it as the subjects variable with the following syntax: MIXED Rt2Adj BY Relatedness Wordness Subjects   /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1)   SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)   /FIXED=Relatedness Wordness Relatedness*Wordness | SSTYPE(3)   /METHOD=REML   /PRINT=DESCRIPTIVES   /RANDOM=Subjects | COVTYPE(VC)   /REPEATED=Relatedness*Wordness | SUBJECT(Subjects) COVTYPE(UN). I get an error when I try to run the analysis ("Iteration was terminated but convergence has not been achieved. The MIXED procedure continues despite this warning. Subsequent results produced are based on the last iteration. Validity of the model fit is uncertain.") I'm assuming that this is because the model is over-specified, because "subjects" was assigned as both the subjects variable and as a random effects variable. So my question is this: in the GUI, SPSS says "Specify subjects for models with correlated random effects". If I specify subjects under this heading, does the model interpret this variable as if it were a random effect variable (even though it doesn't list "subjects" as a random effect under the Model Dimension table the way it does if I manually enter "subjects" as a random effect variable)? Is there any other way to ensure that subjects is entered into the model as a random effect? I am having trouble finding a clear answer to this question because many of the other questions I've seen posted have used a different variable (not subjects) as a random effect. A second question is the following: I used the first syntax pasted above to analyze a very small dataset (4 subjects). Although the analysis ran without errors on a larger dataset, when I ran the analysis on this small dataset I got the following error: "Iteration was terminated but convergence has not been achieved. The MIXED procedure continues despite this warning. Subsequent results produced are based on the last iteration. Validity of the model fit is uncertain." Can a very small sample size cause this error?


title: Estimating good buying and selling prices on eBay using eBay sales history

text: **I am trying to create the algorithm that will tell me the good price to buy stuff and a correct re-selling price based on eBay item history.** The main variables would be the auction type (buy now or bidding), price, if the item was sold and some sort of margin that has to be made on each transaction (in percents or dollar value) This is the spreadsheet with example data: Excel sheet. The first few rows appear below.                AUCTION_TYPE PRICE   DATE_SOLD    SOLD_Y_N      BUY-NOW      $300.00 May-10 22:26 Y      AUCTION      $300.00 May-13 09:11 Y      AUCTION      $300.00 May-13 09:12 Y      BUY-NOW      $474.89 May-14 12:30 N      It is a table of **Auction_Type** (a binary covariate equal to "Auction" or "Buy now"), **Date_Sold** (to the minute), **Price** (consistently in American dollars)--which is an outcome for auctions but is predetermined for "Buy now" sales--and **Sold** , a binary outcome indicating whether the item was sold. **Herein is one of the principal complications:** Price, which is something I would like to estimate, functions sometimes as an "independent" variable and sometimes as a "dependent" variable, depending on the type of auction. This situation therefore does not seem to fit within a standard regression setting. A good solution may require some innovative statistical thinking. Well, I think that I we need to establish some sort of MINIMUM amount of sold/unsold items per week. For example, it is hard to predinct anything if there are 5 records exist. But if we would work with 50 records per week and with the period of coupld of weeks, we could narrow the probability. In addition, autioned items have the higher probability to be sold for lesser value, than the buy-it now. But, on another hand, some unsold auctions can be a warning sign that indicated that this item might not be that valuable. What I am trying to figure out- is there any resonable solution can be created that would allow to buy things based on the item selling/unselling history, and then to be able to resell the same item with the profit and with very high probability of success.


title: Inferential statistics on census data when response rate low

text: I have read the comments about inferential statistics on census data. I am doing PhD and also facing the same problem with a little difference. I have collected census data from the population of deans, directors,chairpersons and head of departments because of lack of availability on their part and can get about 64 % response rate. can i take this 64% response as my sample and by applying inferential statistics eg t test and One way ANOVA draw conclusion regarding mean difference on the basis of gender, age, experience, administrative posts and generalize the results on the whole population. iif it is so then please send me the reference of its justification so that i can defend this conceptual point in my phD thesis. i have been trying a lot but can not find the reference. please is there any one who can help me i shall appreciate and acknowledge it in my dissertation.


title: Test theory standard error

text: Is it true that if you want to make a CI for individuals: S(E)=S(x) and for groups S(E)=S(x)*sqrt(1-rtt)? Thanks in advance!


title: How to correctly assess the correlation between ordinal and a continuous variable?

text: I'd like to estimate the correlation between: An ordinal variable: subjects are asked to rate their preference for 6 types of fruit on a 1-5 scale (ranging from very disgusting to very tasty) On average subjects use only 3 points of the scale. A continuous variable: the same subjects are asked to quickly identify these fruits, which results in an mean accuracy for the 6 fruits. Is Spearman rho the best method to analyze these data and/or are there other good methods I could consider?


title: Application "tuning" based on statistics: Lowest-Learning-Curve API for Java programmer?

text: I want to write an application, where users will have to choose between different "configurations", and I am looking for a way to statistically validate the "balance" of those "configurations", to make the equally "attractive" to users. Ideally, I'd like to use a Java API, as it would make it easy to "round-trip" (the main application will be written in Java) by running the simulation that generates the data, then run an analysis on it, then modify the "configurations", and re-run the simulation, to see how near we are now to the desired state. I've had a solid background in math/statistics, but that was many, many years ago. So I'm looking for something where I can quickly get started, rather then a full-blown, but complex, system, which would require me to spend a lot of time refreshing my math/stats knowledge. To rephrase, I'm looking foremost for an API that allows me to statistically analyse data in Java, rather the "create visualisation" (although, of course, creating visualisation would be useful to quickly validate the analysis code, but the visualisation could be generated by another program). EDIT: Actually, I'm pretty good at Python too, but that would make the round- trip a bit more difficult.


title: Is my understanding of estimation correct?

text: I have two time series, $y$ and $x$, where $y$ is sampled monthly and $x$ is sampled daily. I plan to use a mixed frequency regression(MIDAS) where I assign custom weights to the daily data points in $x$, and by doing so averaging it to $y$'s frequency. My proposed model is as follows: $y = \beta_y L(y)+ \sum_{n=0}^k \beta_nL^n(x) + \epsilon $ where $ \beta_n $ defines a polynomial function, g(n), for some function g. and $L$ is a lag operator such that $L(x_t) = L(x_{t-1})$, $L^2(x_t) = L(x_{t-2})$ One of my supervisors(a professor) for my project commented on this proposal, saying he was " _confused about what you have in mind for my regression_ " He went on to ask me " _There are no parameters in it, for instance, so what would you estimate?_ ". **Will I be right if I told him, that $\beta_y$ and $\beta_n$ are the parameters that are to be estimated, to find out whether we can build a model to forecast monthly $y$ given daily $x$.** Either my concept of estimation is wrong, or the supervisor didnt bother reading through my proposal properly; It has only been email exchanges between us so far, and I havent replied to this one yet. I am confused, please help clarify. Edit: Sorry about $\beta_1$ appearing twice, didnt notice it; I just randomly came up with all those notations, please refer to this answer here for the exact MIDAS regression model ill be using. The polynomial is also explained there. (I was the one who posted the linked question, which is related to the same model talked about in this question)


title: How to measure test set error with logistic regression

text: In measuring the performance of a model, I divide my data into 2 sets, the training set and the test set, fit my model to the training set and then try to predict the results of the test set. If I'm looking at binary classification, I expect to classify my results into 0's and 1's. However, the output of the logistic regression is probabilities. So if I predict a probability of 51%, but I classify this as a 1 because it is > .5, and I get it wrong 49% of the time, isn't my model right as opposed to being 49% wrong? Would it be a better measure to check if my error rate is close to the expected error rate for the model?