title: ARIMA estimation by hand

text: I'm trying to understand how the parameters are estimated in ARIMA modeling/Box Jenkins (BJ). Unfortunately none of the books that I have encountered describes the estimation procedure such as Log-Likelihood estimation procedure in detail. I found the website/teaching material that was very helpful. Following is the equation from the source referenced above. $$ LL(\theta)=-\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(\sigma^2) - \sum\limits_{t=1}^n\frac{e_t^2}{2\sigma^2} $$ I want to learn the ARIMA/BJ estimation by doing it myself. So I used $R$ to write a code for estimating ARMA by hand. Below is what I did in $R$,   * I simulated ARMA (1,1)   * Wrote the above equation as a function   * Used the simulated data and the optim function to estimate AR and MA parameters.   * I also ran the ARIMA in the stats package and compared the ARMA parameters from what I did by hand. **_Below is the comparison:_** ![Comparison](http://i.stack.imgur.com/UUigG.png) **Below are my questions:   * **Why is there a slight difference between the estimated and calculated variables ?**   * **Does ARIMA function in R backcasts or does the estimation procedure differently than what is outlined below in my code?**   * **I have assigned e1 or error at observation 1 as 0, is this correct ?**   * **Also is there a way to estimate confidence bounds of forecasts using the hessian of the optimization ?** Thanks so much for your help as always. Below is the code:               ## Load Packages          library(stats)     library(forecast)          set.seed(456)               ## Simulate Arima     y <- arima.sim(n=250,list(ar=0.3,ma=0.7),mean = 5)     plot(y)          ## Optimize Log Liklihood for ARIMA          n = length(y) ## Count the number of Observatuions     e = rep(1,n) ## Initialize e          logl <- function(mx){            g <- numeric       mx <- matrix(mx,ncol=4)            mu <- mx[,1] ## Constant Term       sigma <- mx[,2]        rho <- mx[,3] ## Ar coeff       theta <- mx[,4] ## Ma coeff            e[1] = 0 ## Since e1 = 0            for (t in (2 : n)){         e[t] = y[t] - mu - rho*y[t-1] + theta*e[t-1]       }            ## Maximize Log Liklihood Function        g1 <-  (-((n)/2)*log(2*pi) - ((n)/2)*log(sigma^2+0.000000001) - (1/2)*(1/(sigma^2+0.000000001))*e%*%e)            ##note: multiplying Log Likelihood by "-1" inorder to maximize in the optimziation       ## This is done becuase Optim function in R can only minimize, "X"ing by -1 we can maximize       ## also "+"ing by 0.000000001 sigma^2 to avoid divisible by 0       g <- -1 * g1            return(g)          }          ## Optimize Log Liklihood     arimopt <- optim(par=c(10,0.6,0.3,0.5), fn=logl, gr = NULL,                      method = c("L-BFGS-B"),control = list(), hessian = T)     arimopt          ############# Output Results###############     ar1_calculated = arimopt$par[3]         ma1_calculated = arimopt$par[4]     sigmasq_calculated = (arimopt$par[2])^2         logl_calculated = arimopt$val     ar1_calculated     ma1_calculated     sigmasq_calculated     logl_calculated          ############# Estimate Using Arima###############     est <- arima(y,order=c(1,0,1))     est


title: In R how to compute the p-value for area under ROC

text: I struggle to find a way to compute the p-value for the area under a receiver operator characteristic (ROC). I have a continuous variable and a diagnostic test result. I want to see if AUROC is statistically significant. I found many packages dealing with ROC curves: pROC, ROCR, caTools, verification, Epi. But even after many hours spent reading the documentation and testing, I couldn't find how. I think I've just missed it.


title: Normality tests for histograms

text: I understand that one quick way to test for normality is to create a histogram from my data and see if it "looks bell-curved". Alternatively, I can use some normality test like Pearson's chi-squared. But what if all the data I have is a histogram? Are there any numerical tests to evaluate if the underlying data is a sample drawing from a normally distributed random variable?


title: Is there more to probability than Bayesianism?

text: As a student in physics, I have experienced the "Why I am a Bayesian" lecture perhaps half a dozen times. It is always the same -- the presenter smugly explains how the Bayesian interpretation is superior to the frequentist interpretation allegedly employed by the masses. They mention Bayes rule, marginalization, priors and posteriors. What is the real story? Is there a legitimate domain of applicability for frequentist statistics? (Surely in sampling or rolling a die many times it must apply?) Are there useful probabilistic philosophies beyond "bayesian" and "frequentist"?


title: Any bibliographical reference where it says that the Dice Coefficient is the same as the F1 - measure?

text: Running tests on a segmentation algorithm I realized that the Dice coefficient was consistently giving me the same results as the F1-measure. I want to mention this in a paper I'm writing, that I've only found this as reference (and I think it looks bad to put into a paper). Any other sources?


title: Problem simulating AR(2) process

text: **OP EDIT** : There where no problem with this. The problem was with the method I was using for obtaining the PACF. Apparently it doesn't work quite well in this case (I was using the scikits/tsa python package to obtain the PACF via the YW equations). Testing the coefficients in R worked like a charm. I'm trying to simulate an AR(2) process, but I seem to be getting awful results. The way I'm doing it is as follows: if I want to simulate 1000 points of an AR(2) process with coefficients, say, a1=0.1 and a2=0.5, I simulate a realization of 2000 points of a white noise process (in my case I simulated 2000 points drawn from a normal, ~N(0,1) distribution), where I'll use the first 1000 points as burn-in points. Suppose I store this realization in a vector W[t]. Then, I simulate the actual AR(2) process by iterating on a new vector, X[t], as follows:               X[0]=W[0]          X[1]=W[1]          X[2]=a1*X[1]+a2*X[0]+W[2]          X[3]=a1*X[2]+a2*X[1]+W[3]          ...          X[i]=a1*X[i-1]+a2*X[i-2]+W[i]      Finally, I return the last 1000 values of the vector X[i]. The problem is that when I plot the actual Partial Autocorrelation Function (i.e. when I estimate the coefficients of the AR(2) process I generated), I get wrong coefficients for a1 (I get an acceptable value for coefficient a2, though). What am I missing? Here's my Python function for the simulation:               from numpy.random import normal     from pylab import *          # INPUT:      # a:        Is the array with coefficients, e.g. a=array([a1,a2]).     # sigma:    The white noise (zero-mean normal in this case) standard deviation.     # n:        Number of points to generate.          def ARgenerator(a,sigma,n,burnin=0):            if(burnin==0):         burnin=100*len(a) # Burn-in elements!       w=normal(0,sigma,n+burnin)       AR=array([])       s=0.0       warning=0       for i in range(n+burnin):           if(i<len(a)):             AR=append(AR,w[i])           else:             s=0.0             for j in range(len(a)):                 s=s+a[j]*AR[i-j-1]             AR=append(AR,s+w[i])       print 'Measured standard deviation: '+str(sqrt(var(w[burnin:])))       return AR[burnin:]      **Edit:** by an MCMC simulation (where I use the autocovariance matrix for an AR(2) process with a multi-variate gaussian likelihood), I get right the value of sigma (the std. deviation of the white noise process) and the a2 coefficient. However, the a1 value I obtain has nothing to do with the true one.


title: What's the statistical way to rainbow table a web service?

text: Rainbow tabling a web service is getting a whole bunch of pages by guessing their IDs:               /item?ID=1     /item?ID=2     ...     /item?ID=n      Ideally, I want to intelligently detect the max ID and set `n` to that. Is there a statistical approach here? I was imagining like a "comb" where I check a radius of IDs around an exponentially increasing value:               9,10,11     99,100,101     999,1000,1001     9999,10000,10001      or maybe the exponentiation is dangerous, so just step by a constant:               999,1000,1001     1999,2000,2001     2999,3000,3001      E.g. if none of `2999,3000,3001` are valid IDs, then I set `n` to `3000`. Should I be using some crazy distribution or formula with _e_? _Note:_ when DB admins want to make their IDs "unguessable," they often just choose a big number to start them at - say 100000, so item 1 is 100001, item 2 = 100002, etc. _Note:_ the brute force approach is just to set n to `2147483647` (the max value of MySQL's `INTEGER` type) or a smaller value if you know roughly how many objects there are in the database, and go one-by-one.


title: Combining errors

text: I've got this tricky problem to solve. Hope somebody can help me. Let's suppose I want to know the mean radius of the trees in a wood. For each tree i I measure the radius at the top, at tha center and at the bottom, and i get 3 values with the systematic uncertain given by the resolution of my instrument (always the same): $X^i_{bot}$ $X^i_{mid}$ $X^i_{top}$ $\Delta X_{sys}$ How do I combine the values and how o I get the errors? My first approach was to combine the 3 measures statisticaly in order to get a mean $X^i$ and a st. deviation $\sigma$, and then add the systematic error. $\Delta X^i = \sqrt{\sigma^2+\Delta X_{sys}}$ Then to combine the different trees I've made a mean of central values and combined the error with the formula $\Delta X=\sqrt{\Sigma (\frac{\delta f}{\delta X_i} \Delta X_i)^2}$ But I don't think that wuold be correct. Furtheremore, I don't only need the trees mean radius, but I need the mean radius of a ring of a tree of that forest, provided that each tree has a different number of rings $n^i$. I really don't know how to make this last part.


title: Calculating the threshold of a set of values

text: **Question:** The context of this question is actually finance, however the question itself is a statistical issue. Suppose I have the following expression: $$\rho = \frac{2\bar{x}}{(s^*_x)^2}+1 .......... (1)$$ where $\bar{x} = \frac{1}{T} \sum_{t=1}^{T} x_t$ $(s^*_x)^2 = \frac{1}{T} \sum_{t=1}^T (x_t - \bar{x})^2$ Assume $T$ is some fixed constant. Note: $T$ is just the total number of observations of $x_t$. I have data on $x_t$ for each 'entity' (here an entity just simply refers to a firm/company). In total, I have 2228 entities and for each entity I have $T$ observations of $x_t$. For each entity, I substitute the $T$ observations of $x_t$ into Eqn. $(1)$ and obtain a value for $\rho$. Thus in total, I have 2228 values of $\rho$. Now, a large value of $\rho$ means the entity is "bad" and a small value of $\rho$ means the entity is "good". However, the problem is how large does a value of $\rho$ have to be in order to classify an entity as "bad"? That is, what is the threshold such that if $\rho$ exceeds the threshold value, then we can classify the entity as "bad"? For example, let's say the threshold is $400$, if entity A had a $\rho = 300$ while entity B had a $\rho = 1000$, then entity A is "good" while entity B is "bad". * * * **My attempts so far:** My first attempt was to try get data on an entity that is known to be "bad" and then calculate its $\rho$ and use this value as the threshold. The problem is that I cannot obtain data on "bad" entities. For my next attempt, I obtained the empirical distribution by applying a kernel density estimator on the 2228 values of $\rho$. Then I calculate the 99th percentile (for robustness, I also calculated the 97.5th and 95th percentile) of this pooled distribution and use this value as the threshold. However, the main critique is that this is too arbitrary and there is not enough rationale for using this method. * * * **Queries:** So I am wondering if anyone has any ideas on how what statistical/mathematical techniques/methods I can apply to derive appropriate thresholds for $\rho$. Currently, I really have no idea on what tools are available for this problem.


title: Which machine learning algorithms can be scaled using hadoop/map-reduce

text: Scalable machine learning algorithms seem like the buzz these days. Every company is handling nothing short of _big data_. Is there a textbook which discusses what machine learning algorithms can be scaled using parallel architectures like Map-Reduce, and which algorithms cannot? Or some relevant papers?


title: Need some advice on data analysis

text: I collected data from 500 participants. Here's an example of the data: Response 1:   Age: 25   Gender: Male   Choice 1: 1   Choice 2: 1   Choice 3: 0   Choice 4: 1   Choice 5: 1   ... Up to choice 30 Response 2:   Age: 20   Gender: Female   Choice 1: 1   Choice 2: 1   Choice 3: 0   Choice 4: 1   Choice 5: 1   ... All questions were categorical and the choices were binary outcomes. My goal:   a) To see what choices were preferred in general (using the mean, since the data is normal)   b) Does gender make a difference in choice (using linear regression) What I really want advice on:   1) I put in control variables and it shows a preference for choice 1 by about 3%. What should I do to correct this to make it 50/50? Or should I even do that at all?   2) There are responses to the survey where participants just selected 1 for all responses. Should I reject these? Thank you!


title: Confidence Interval for a Random Sample Selected from Gamma Distribution

text: Working on a homework question and having some trouble... Any help would be greatly appreciated. Based on a sample 1.23, 0.36, 2.13, 0.91, 0.16, 0.12 from the GAM$(2,\theta)$ distribution, find an exact 95% CI for parameter $\theta$. So we know GAM$(\alpha, \lambda)$ has the pdf $f(x)= \dfrac{\lambda^{\alpha}}{\Gamma{(\alpha)}} x^{\alpha - 1} \ e^{-\lambda x} $. Therefore our random sample is distributed with pdf $f(x)=\theta^{2} x e^{-\theta x}$. I understand that because the question asks for an _"exact"_ confidence interval, that I need to find the pivotal variable. The problem I am having is that most examples I find are along the lines of a random sample... $X_1,...,X_n \sim N(\theta, \sigma^{2})$ if $\sigma$ is known then $Z= \dfrac{\bar{X}-\theta}{\frac{\sigma}{\sqrt{n}}}\sim N(0,1)$, is pivotal. And from there finding the CI is relatively simple. I guess I am at a loss as to how one would go about finding the pivotal variable when things are not normally distributed. Thank you for your help, any suggestions would be appreciated.


title: Bias in predicting survival with time varying covariates

text: I'm trying to predict survival probabilities with time-varying covariates. My dataset constitutes a variety of subjects who enter the study at different dates and receive multiple follow-ups. For each subject, I count time from the date of entry into the study (rather than from the beginning of the study). When I evaluate my predictions on the test set, I noticed that the model performs especially badly for people who enter the study shortly before it ends. Is there something in the way that I have designed this model that is introducing bias, or could this simply be because I have less information on later entrants? * * * In case this is helpful, I'm using the R survival package to fit a cox model with `Surv(start,stop,event)` form, where I have multiple observations for each individual. I'm using the predictSurvProb in PEC to produce the predicted survival probabilities for my test set.


title: Confusion related to a post related to statistics in xkcd

text: I came across this post http://xkcd.com/1132/ and this image ![enter image description here](http://i.stack.imgur.com/g8yXV.png) I didn't get why the Bayesian statistician is saying no. Any suggestion


title: Integral of a conditional uniform distribution leads to improper integral

text: I have two uniforms distributions, $X_1 \sim\it{U}(a,b)$ and $X_2\sim\it{U}(X_1+\delta,b+\delta)$. I would like to compute $P(X_2\in[a+\delta,b+\delta])$. So I do this: $$\begin{eqnarray*} P(X_2\in[a+\delta,b+\delta]) & = & \int_{a}^{b} P(X_2\in[a+\delta,b+\delta]|X_1=y)\cdot P(X_1=y) dy\\\ & = & \frac{1}{b-a} \int_{a}^{b} \frac{1}{b-y} dy\\\ & = & \frac{1}{b-a} \ln\Bigg(\frac{b-a}{b-b}\Bigg) \end{eqnarray*}$$ The wrong result is related to the fact that a uniform distribution has probability $\frac{1}{x}$, but it doesn't work for $x<1$ (in this case, the fraction becomes bigger than 1...). How do we solve this problem formally, in a clean way?


title: Variables importance: who can do the most pushups?

text: I don't know enough math to formulate an intelligent question on this so I'll give an example. I'd like an answer to my example but also I'd like to know the jargon I need to be able to research it further. Suppose you have a large set of data of information about people: their height, weight, age, and number of pushups they can do. How would you go about finding:   1. Which height, weight, and age combinations can do the most pushups?   2. Given a height, weight, and age, what is the probability they can do x number of pushups? Ultimately I'd like to know, given a set of 10 people's heights, weights, and ages, who would statistically be the most likely to do the most pushups.


title: How many samples account for all the variance?

text: I have a data set of 100,000. If I want to pick a subset (sample), how can I ensure that it's a "good" sample? I know that I can make it unbiased by picking at random, but how do I know how many to pick to account for all the variance in the dataset?


title: How is it possible that this ratio revert to the mean?

text: I'm testing the ratio on the chart below with ADF. Ratio chart The function I'm using is:               adfTest(ratio, type='nc')      on fUnitRoots library. The result (p-value) is: **0.01301** It is **very low** so the serie should be mean-reverting. (@Statistical Please do not hate me for this brutal explanation) Seeing the chart the ratio doesn't seem very "mean-reverting", I'm talking about a visual check of the ratio. What do you think about that? Thank you!


title: Effect of sample size in t test?

text: > **Possible Duplicate:**   >  Why does t statistic increase with the sample size? In t test,the test statistic is $$ t = \frac{\overline{x} - \mu_0}{s/\sqrt{n}} $$ I wonder what effect of the sample size $n$ has on the t test? For example, as $n$ increases, I thought $t$ will increases at first, but later I realized $s$ is at the scale of $1/\sqrt{n-1}$. So $s/\sqrt{n}$ seems not to decrease as $n$ increase? What other effect does $n$ have? Thanks!


title: Value of alpha and beta in Holt's exponential smoothing method

text: How to choose the best values of alpha and beta in Holt's exponential smoothing? Leaving it upon R gives me $\alpha$ =1. Is this appropriate? Entering different values of alpha and then comparing with the real data shows best result for $\alpha$ = 0.45. But then R calculates $\beta$=0.99. Is this fine?


title: Alternative to block bootstrap for multivariate time series

text: I currently use the following process for bootstrapping a multivariate time series in R:   1. Determine block sizes -- run the function `b.star` in the `np` package which produces a block size for each series   2. Select maximum block size   3. Run `tsboot` on any series using the selected block size   4. Use index from bootstrap output to reconstruct multivariate time series Someone suggested using the meboot package as an alternative to the block bootstrap but since I am not using the entire data set to select a block size, I am unsure of how to preserve correlations between series if I were to use the index created by running `meboot` on one series. If anyone has experience with meboot in a multivariate setting, I would greatly appreciate advice on the process.


title: Data mining of time series

text: I have a dataset which consist of n time series variables $X_1$..$X_n$ , and a time serie output $Y$. I would like to mine the data to find if some functions (lagged or not) of the $X_i$ can predict $Y$. The $X_i$ may have missing values (the data is unbalanced). Which algorithms would you suggest to do this?


title: How to calculate  95% confidence intervals for  partial eta squared and eta squared?

text: Following on from this question on reporting eta squared, I would like to report 95% CIs for my effect size point estimates (partial eta and eta squared values).   * What statistics do I need to know and is there a straight forward formula for calculation?


title: Getting started with multivariate testing

text: I'm looking for the best combination of variables to result in a click to the next page of my website... For example the size and color of a submit button. I have a system that currently tests each variable independently to find the best value for that variable. I would like to take the next step to testing combinations of these variables. I'm looking for help with 2 steps in this process:   1. How to save the result sets - Currently I store each variable independently in mysql with no link to the other variable values that were tested. What are some best practices for recording this data (schema design)?   2. How to analyze the result sets - Stats is not my strong suit, and I don't have any previous experience with stats software... I'd like a more-or-less black-box type of approach.... I supply the raw data (from step 1) and the analysis software spits back out to me the best combination of optimized variables.


title: Why don't my two-stage least squares results make sense?

text: I used SPSS 19's 2SLS procedure (which is very straightforward, with almost no optional specifications) to predict Y from X after X was predicted based on I, an instrumental variable. Then I tried to match those results by running 2 separate OLS regressions. First I obtained the predicted X values by regressing X on I. Then I regressed Y on this "predicted X" variable. The results hardly matched those of the integrated 2SLS procedure. The _b_ matched, but not RSQ or Beta (which was wildly different) or _t_. Why would this happen? I made sure to filter by the same subset of cases for each procedure.


title: Basics of extreme values / high-water marks?

text: With real-valued $X_1, X_2, \ldots$, define   $Max_n := \max(X_1,\ldots,X_n)$ record value or high-water mark   $NextMax_n :=$ the next greater high water, $Max_{n+m} > Max_n$   $Up_n := NextMax_n - Max_n$. Can anyone suggest an introduction to estimating $NextMax$ and $Up$, for   1. $X$ i.i.d. from a known distribution,    2. nonparametric, estimate-as-you-go?  Either a book chapter or an online stats course would do nicely.


title: Sphericity assumption in multivariate approach

text: I read a few texts about the multivariate approach for repeated measures where k repeated measures are transformed to (k-1) "difference scores". All texts say that using this approach, the assumption of sphericity is not required. But none of those texts explains why this is. Can anyone explain why? Thanks


title: Showing that $s^2$ is an unbiased estimator of $\sigma^2$

text: $s^2 = \sum \frac{(x_i - \bar{x})^2}{n-1} $ which apparently equals $ \frac{\sum(x_i^2) +n\bar{x}^2 - 2n\bar{x}^2}{n-1}$. Does this just come from expanding the numerator and using the fact that $\bar{x}$ (the average) is not dependent on $i$? How might one continue this to show for certain $E(s^2) = \sigma^2$?


title: Square distance and likelihood in k-means

text: In k-means algorithm, the distance minimization step is equivalent to maximize likelihood: $P(X|\theta)$ or to maximize posterior distribution $P(\theta|X)$?   I think it's more logical to maximize posterior distribution $P(\theta|X)$ because it corresponds to classify $X$ in a cluster and usually it is posterior distribution $P(\theta|X)$ which is used for that. If it corresponds to maximize likelihood $P(X|\theta)$ can you explain me why?   Thank you


title: Merging different datasets for binary classification

text: I have **four** ($D_1, ... D_4$) different datasets from which I want to learn a robust _binary_ classifier. I am currently using `SVC` and `LinearSVC` from the `sklearn` library. The problem is that only one of this dataset is already split in training and test test. The others three are not. Albeit the feature set is essentially the same for all of them, each dataset has a different number of samples and they appear to be non coherent with respect to _"quality"_. To clarify in dataset $D_1, D_2$ there are some positive samples that are incorrectly labeled with class $-1$. What is the best way to merge the datasets into a single one trying not to introduce any bias?


title: Performance comparsion study on real data in MATLAB (machine learning)

text: I need to compare performance using 2-class classifiers–an LDA classifier and Exact Bayes in MATLAB. I have to use this dataset. Can anyone give me any advice how to do that (at least the steps of the process)?


title: Sum of squares of residuals instead of sum of residuals

text: As I was working my way through my book about statistics, I came across the topic of linear regression. During the chapter the author begins with explaining that you want to minimize the residuals in order to make your y = a + bx as good a fit as possible: I do understand this, but halfway the chapter all of a sudden the residuals change into sum squares of residuals. Why is this done? I have been googling, but couldn't find the right answer. Who would like to help me understand why the sums of squares of the residuals is used instead of just the sums of residuals? Kind regards, Bas


title: Why we need variance in this world?

text: Why we need variance in this world? What is the purpose to make such function and how does it work? I know that it is a measure of how the data spread, but why we don't just use the absolute variation to solve this problem?


title: Squared Error vs Absolute Error loss functions

text: The two most popular types of loss functions are 1) squared error: $L(y,f(x))=(y-f(x))^2$ --> best estimate is the $E(Y|x) $   2) absolute error: $L(y,f(x))=|y-f(x)|$ --> best estimate is the $median(Y|x)$ I have two questions. 1) Why do people use the squared error method? The absolute error method makes much more intuitive sense. You get the difference between the actual and the estimate. Plain and simple. If you square the difference, then won't you get "warped" values depending on the size of the difference? 2) This also got me thinking about what is "expected value." Expected value is defined as the mean. However, the best estimate under the absolute error loss function is the median. So is the "expected value" the median? I would very much appreciate it if someone can help me clarify my thinking. Thanks.


Why are the squares of the vertical distances taken rather than the absolute values?


title: Why squared residuals instead of absolute residuals in OLS estimation?

text: Why are we using the squared residuals instead of the absolute residuals in OLS estimation? My idea was that we use the square of the error values, so that residuals below the fitted line (which are then negative), would still have to be able to be added up to the positive errors. Otherwise, we could have an error of 0 simply because a huge positive error could cancel with a huge negative error. So why do we square it, instead of just taking the absolute value? Is that because of the extra penalty for higher errors (instead of 2 being 2 times the error of 1, it is 4 times the error of 1 when we square it).


title: Why is the variance squared, and does it mean the same thing as the standard deviation?

text: I am still confused, despite similar questions being asked, about the difference between the variance and standard deviation in statistics. Why is the variance squared?


title: Why square the difference instead of taking the absolute value in standard deviation?

text: In the definition of standard deviation, why do we have to **square** the difference from the mean to get the mean (E) and take the **square root back** at the end? Can't we just simply take **the absolute value** of the difference instead and get the expected value (mean) of those, and wouldn't that also show the variation of the data? The number is going to be different from square method (the absolute-value method will be smaller), but it should still show the spread of data. Anybody know why we take this square approach as a standard? The definition of standard deviation: $\sigma = \sqrt{E\left[\left(X - \mu\right)^2\right]}.$ Can't we just take the absolute value instead and still be a good measurement? $\sigma = E\left[|X - \mu|\right]$


title: Using results of regression on raw observation values to approximate results of regression on relative change between observations (Simple, Linear)

text: this is my first time on Stack Exchange so if I did something wrong please tell me. I have a time series dataset.   There is an observation $(y,x)$ for each continuous time $t$.   Let’s say for each day since 2014-01-01, i.e. $t \in \left\\{0, 1, 2, …, 108 \right\\}$, where $108$ means today. I ran a regression on the model   $y = \alpha + \beta x$   and found $\alpha$ and $\beta$ with $R^2 > 0.975$ My question is, can I **approximate** $\beta'$ in the model   $y' = \alpha' + \beta'x'$   where   $y'=\frac{y_{t+1}}{y_t}-1$   $x'=\frac{x_{t+1}}{x_t}-1$   **without** running regression on that model, just by using $\hat{\beta'}= \frac{\beta \times x_{108}}{\hat{y_{108}}}$ where $\hat{y_{108}}$ is estimated using the original model. I'm curious because I found that my $\hat{\beta'}$ is usually within 1% of $\beta'$ Even a simple Yes or No with help me a lot. Thanks!


title: Do skewness and kurtosis uniquely determine type of distribution?

text: Inspired by this answer, I have following question: Is it enough to know just skewness and kurtosis in order to determine distribution that data comes from? Is there any theorem that implies this? Moreover, can we replace skewness and kurtosis by any other pair of moments (for example expected value and variance)?


title: Cronbach's Alpha for open questions

text: In the questionnaire I asked respondents from two countries how many job offers they received from 5 sources in the last 6 months. There are 5 questions - one for each source. It is an open question, without a scale as the two countries strongly differ in labor market conditions. The range of responses is from 0 to 30. Data is very skewed (most people indicating 0). Can I calculate Cronbach's Alpha for each country separately? Values are very low then. The values are much better if I transform the data into a scale (0, 2-3 offers, 4-5 offers, 6-9 offers, more than 10 offers). What is the best approach?


title: Inflation as an independent variable

text: Assume a model like this, basically explaining stock market returns with a bunch of stuff:               stockReturn(t) ~ bondReturn(t) + moneyMarketReturn(t) + inflation(t) + somethingElse(t)      1) Does using inflation as an independent variable bring any significant problems? 2) Should the results differ from using real returns instead, in this case using the inflation to calculate real stock, bond, and money market returns (assuming that the somethingElse is a non-economic variable). What would affect the difference? So here we would do this for the stock, bond, and money market returns to get the real returns:               realReturn(t) = ( 1 + return(t) )/( 1 + inflation(t) ) - 1      and then drop inflation from the model.


title: How to examine differences between two groups across three different experimental conditions?

text: I have data for two groups. The 2 groups are matched on a number of variables, but differ in the scores on an alcohol measure (one group is high, and one is low). However the groups have been created so that they capture the extremes of the alcohol measure (all very low in one group, all very high in another) therefore the alcohol measure that separates them is not strictly continuous. I then have data for the groups on task performance in 3 different conditions: reward, neutral, and punishment. I want to compare the groups for their task performance in each of these conditions, and additionally see if each group responds differently across the conditions. I have conducted t-tests investigating if the means are different for each group, in each different condition. The groups significantly differ in one of the conditions. Question 1: Would it be valid to do a RM-ANOVA treating each group as independent of each other, thus using a split file in SPSS and comparing scores across the three conditions within each group? Question 2: I am considering a MANOVA, using task performance in the 3 conditions as dependent variables and group as the independent variable. If I were to do this, would it replace the t-tests or can it be used as additional examination (whereby I would have data using a univariate approach in the t-tests, and then using a multivariate approach using the MANOVA). Any thoughts?


title: Is it appropriate to report means and standard deviations for composite scores?

text: **Following the posthere, I am wondering if reporting the mean and standard deviation is appropriate when describing a composite score.** My composite score is the total of responses on several likert scale items and count variables. I have made sure that all are on a five point scale, based on their level of importance, so as to preserve variability. This is the weighting system I have used to assign the scores:   * 0 = of no importance   * 1 = of some importance   * 2 = of medium importance   * 3 = of high importance   * 4 = of absolute importance  I have then added all scores together to arrive at the composite score. As such, this is based on a formative measurement model framework. I am aware of issues relating to the purity of scales and I am hoping my logic of using "level of importance" as the weighing criterion would provide the necessary meaning I need to interpret the composite score. For example, if I get a composite score of 3, I know it is of high importance. (I apologise for my non-statistical language.)


title: How to pick a good threshold for multimodal histograms

text: I am trying to do multi-modal detection in a vote map (histogram) however, side lobes obviously cause some false peaks. So using mean shift will not do as it will find all of the local maxima. So I thought I could first threshold the histogram and then apply mean shift. However, I am struggling to find a suitable threshold value (derivable rather than a hard value). This is because the data can range from very close to zero to about 19 or 20. It is a 480x640 vote space. For this reason, the mean is inappropriate as it is usually a very small number as the number of tiny peaks is generally much higher than the true peaks. Similarly standard deviation is quite small. Any suggestions are welcome. Thanks


title: Formula to generate Critical t-values for t-test (instead of using a look-up array)

text: I am searching for the formula that will generate the critical values in a t-test instead of using a static critical value table. I am programming in Java and not R so the qt() function is not what I am looking for. Any help is appreciated.


title: Analyzing the added effect of an individual variable having fitted one variable already, using a generalized linear model

text: I am using a generalized linear model to analyse my data, there are 6 explanatory variables. I have added all different variable combinations in different models, and ranked them according to their AIC weights. Only one variable is really significant in any of the models, but different models with different combinations of variables rank reasonably closely behind the 'best' model with the highest AIC weight which only contains the explanatory variable that is significant. I would now like to analyse the effect of the individual variables a bit further, in the sense that I would like to analyse what the effect is of the non-significant variables, while controlling for the significant one. So given the effect that variable $X_1$ (the significant one) has on $Y$, how much more effect does $X_2$ have for example? I initially thought of fitting a generalized linear model to my data, using only $X_1$, and then fitting another generalized linear model using $X_2$, on the residuals from the first model, however, that method has some major disadvantages I understand. Does anybody have some advice on how else I can do this?


title: What does intercept actually mean in a simple linear model?

text: Given y = a + b*x, does "a" represent simply an intercept or mean y? For y = a + b*ln(x), is there a need to interpret the intercept "a?"


title: R: Problem with NaNs in gamma regression using glm

text: I have found a problem with `NaN`s when I try to fit the model from example 5.1 in the book _Generalized Linear Models. With Applications to Engineering and Life Sciences_. The authors fit a `glm` with gamma distribution and a log link to a dataset about resisitivity (y) and 4 predictors (x1..x4). First, they use a full model: `y ~ x1*x2*x3*x4` and then remove the non significant interactions to arrive to the final model: `y ~ x1+x2 +x3 +x4 + x1*x3 + x2*x3 + x3*x4`. When I try to replicate the full model with R (they use SAS) I get a warning: `dgamma(y, 1/disp, scale = mu * disp, log = TRUE) NaNs produced` and in the `summary(model)` the `st.err` and next ones columns are with `NA`s. Therefore, the significance of coefficients is lost.               x1 <- c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)     x2 <- c(-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1)     x3 <- c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1)     x4 <- c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1)     y  <- c(193.4,247.6,168.2,205,303.4,339.9,226.3,208.3,             220,256.4,165.7,203.5,285,268,169.1,208.5)          res <- data.frame(x1, x2, x3, x4, y)     mod <- glm(y ~ x1*x2*x3*x4, family=Gamma(link="log"), data=res)     summary(mod)          Deviance Residuals:       [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0          Coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept)  5.26476         NA      NA       NA     x1           0.24705         NA      NA       NA     .     .     .     x1:x2:x3:x4  0.36586         NA      NA       NA          (Dispersion parameter for Gamma family taken to be NaN)      Is there some way to circumvent this problem? When I fit the final model: `y ~ x1+x2+x3+x4 + x1*x3 + x2*x3 + x3*x4` the results are exactly the same than in the book. So, why SAS can resolve the full model and R returns errors?


title: Normalized Variation of Information (NVI) visualization of distances to several references

text: I am using NVI (Normalized Variation of Information) to measure the quality of clustering obtained ( a family of them) compared to the Ground Truth and also to the result obtained with another method. Since NVI is a metric, it satisfies the triangular inequality and hence I thought it would be nice to plot the result graphically. For that, I need a tool that given the distances among points can plot them in such a way that all the distances are respected among the different reference points. The graph should be 2D. Do you know of any tool that can do like that? i.e. receiving a matrix of distances The matrix will not be complete, that is, the distance between all the partitions is not calculated. All the distances to the 3 references (2 GT and the additional method) and the distances among the 3 references are calculated. If not a tool, using Matlab or Python and doing some programming is also fine for me. I have around 200 points to plot. EDIT: I have found this answer in SO: http://stackoverflow.com/questions/13513455/drawing-a-graph-or-a-network-from- a-distance-matrix but it doesn't work properly, even for 7 nodes (that can be visually verified the correctness")


title: multiple linear regression with normalization - how to get non-scaled full covariance matrix

text: I am doing a quite complicated multiple regression modelling in physics and have a problem **how to got back to covariance matrix for non-normalized parameters**. I don't know how to calculate the error for the intercept and get the full covariance matrix (that is a matrix that will include a row and column also for intercept). I have a model: $$ d=Gm $$ where $d$ is a vector of my experimental data, $G$ a design matrix and $m$ a vector of model parameters. I need to use a TSVD for regularization (problem is ill-posed) and hence need a normalization. This is done by: $$ d_s=(d-\bar{d})/\sigma(d) $$ $$ Z[i]=(G[i]-\bar{G[i]})/\sigma(G[i]) $$ where G[i] is a i-th column of G. As I removed the mean from $d$ I don't need intercept in scaled model. Therefore the solution of my problem is: $$ m_s=(Z^TZ)^{-1}Z^Td_s $$ where the vector $m_s$ do not include intercept. I can now calculate the covariance matrix for $m_s$ and got back to $m$. I know how to calculate the intercept. This I do with: $$ m_0=\bar{d}-\sum\bar{G[i]}m_s[i] $$ however I have no idea how to get the error for it and calculate the full covariance matrix for $m$. I need this as I need to see it in true physical units. Try to find it by comparing with the results of regression without normalization but couldn't guess anything. Also hard to find in regression books. **Any help will be very appreciated!**


title: How can stochastic gradient descent avoid the problem of a local minimum?

text: I know that stochastic gradient descent has random behavior, but I don't know why.   Is there any explanation about this?


title: Inverting non positive definite covariance matrix

text: I have an expression for a covariance matrix $C$ in terms of the indices $i$ and $j$. In this way I can analytically calculate the elements of my covariance matrix, however when I try to invert $C$ matlab gives a warning about the matrix being close to singular. The inversion therefore doesn't work, by which I mean the multiplying $C$ by the inverse given i do not get the identity. I have tried calculating the pseudo inverse but this also does not work. I have also tried adding a small constant along the diagonal but again the results do not work. In general I'm working with matrices of dimension 1200 but I will give low dimensional example matrix that has the same properties, i.e. the matrices are symmetric about the diagonal and the anti-diagonal                  19.9939  19.9954  19.9958  19.9951  19.9933  19.9905  19.9865        19.9954  19.9973  19.9981  19.9978  19.9965  19.9940  19.9905        19.9958  19.9981  19.9993  19.9995  19.9985  19.9965  19.9933        19.9951  19.9978  19.9995  20.0000  19.9995  19.9978  19.9951        19.9933  19.9965  19.9985  19.9995  19.9993  19.9981  19.9958        19.9905  19.9940  19.9965  19.9978  19.9981  19.9973  19.9954        19.9865  19.9905  19.9933  19.9951  19.9958  19.9954  19.9939      As mentioned in the title the matrix isn't positive, however the the negative eigenvalues are very small suggesting that the matrix is not positive only due to machine precision. The negative eigenvalues are $-1.4048e-14$ and $-2.4571e-15$. How can I go about inverting these matrices?


title: Calculate the arithmetic average and standard deviation of a survey response

text: If I have a market research with population of 100 individuals:   1. 10 spend less than 10 dollars in clothing    2. 20 spend between 10 - 20 dollars in clothing    3. 30 spend between 20 - 100 dollars in clothing    4. the remaining do not spend at clothing.  _How do I calculate to arithmetic average and standard deviation of the data?_


title: Compute mean from discrete normally distributed data

text: I am looking for an estimator of the mean of normally distributed data (with known variance) in the regime where the sampling grid is much coarser than the variance. I.e. in the worst case consider count data where only two bins contain data points. For example assume I sampled 40 data points from a normal distribution with unit variance and have two bins: 1) Bin one from -1.5 up to 0.5 contains 28 data points 2) Bin two from 0.5 up to 2.5 contains 12 data points


title: Accuracy of computing moments from an histogram

text: Let's say I have an histogram of some continuous data with unknown and arbitrary pdf. I don't have access to the original raw data. So I don't know that 3 of my data points are 21, 25 and 27, all I know is that 3 of my data points are in the bin 20-30. How accurate is it to compute the moments of this pdf from the binned data? Can I estimate the error I'm making?


title: Markov chain getting stuck due to insufficient data samples

text: There is a lot of theory on Markov models and output generation, but I cannot locate any information on models getting stuck. I'm trying to create a model of a data set using a Markov model. The data can look like this "abc abb acc baa bcc...", and I want to make an n-gram model. Accordingly, I sample the data set at random, so I get a model like this (example of 2-gram model):   * abc abb -> acc with probability p1   * acc baa -> bcc with probability p2   * ... The problem occurs, when I try to generate an output from the model. Say I initiate the model like this:   * First: abc abb => acc, so the output is now "abc abb acc"   * Second (taking the last two words of the output): abb acc => ??? The model gets stuck, because the data set is not complete, and therefore does not cover every possible combination. When making the mode, the sample "abb acc" was never reached, and thus the output cannot be determined. I initially asked the question on the C.S. forum site (link), and was advised to ask it here. From the C.S. site, smoothing functions like the Laplace smoothing and Good-Turing smoothing were suggested. Is that the best way to go, or does other methods exist?


title: Statistical Design for Categorical Data Analysis with Selection on DV

text: We are interested whether experimental units belong to class A given that they are members of class B; clearly the classes are not exclusive. We have many millions of experimental units that we can sample from at will. The problem is that we only have only discovered 50 examples of class A; there are likely others in our database, but they were not discovered at the time, and now it is impossible to test whether or not they are members of class A. We have every reason to believe that Class A is very, very small, though, so it is unlikely that any given shipment that we have identified as Not-A is, in fact, A. The ordinary research design would simply be to draw a random sample and place them into a 2x2 table according to whether they are A and B, A and Not-B, etc. In this case, however, the probability that a sample will contain any of the examples we have identified as A is vanishingly small, and even if collect a very large sample and we manage to sample a known example of Class A, we not get great results from the chi-square test as the smallest expected cell count will be much closer to 0 than 1; Fisher test requirements of fixed marginal totals are not met. Instead, I propose sampling from the Not-A population and recording B. The resulting table will tell us Pr(B|A), and using our known population of A, and its information about B, to construct the 2x2 table. Using Bayes' rule, we know $Pr(A|B)=Pr(B|A)*Pr(A)/Pr(B)$. Pr(A) could be estimated as a fraction: the number of example of class A in the database over the number of examples. Estimating the number of examples of class A, discovered and undiscovered, can be achieved by hierarchical Bayesian techniques, at the cost of additional assumptions on my part. $Pr(B)$ is trivial for me to estimate, as we know that information for every example. I have several questions: 1\. Is this reasoning sound? 2\. Am I re-inventing the wheel? Are there other, more established methods to deal with situations like these? Any insight or feedback is much appreciated.


title: Nonparametric or robust alternative to factorial repeated measures design

text: **I have a repeated measures design with two factors** (two levels each). Each subject has a number of trials with each of the four combinations of the two factors. Though some combinations weren't carried out on a few of the subjects. **The response variable is proportion data** with lots of zeros, meaning a general linear mixed model/repeated measures anova is not appropriate (as residuals are non-normal). I have tried a generalised linear mixed model, using the _lmer_ function in R ( _lme4_ package), with my two factors as fixed variables and subject as a random variable, but the residuals are not particularly normal. This is true using a poisson error distribution with a log link or using a binomial error distribution with a logit or probit link (though the latter is best). **I was wondering if I could use a robust or non-parametric test in R to look at the effects of my factors on the reponse variable (taking into account the repeated measures)?** I came across Brunner's _nparLD_ package, but it doen't work for my data because a) _each_ subject experiences _each_ combination of factors, and b) there is no time variable.


title: What's a good Nagelkerke value?

text: My model has Nagelkerke $R^2 = 0.282$ as outputted by SPSS. Having read some literature, it claims that this should be counted as a good value. However, since the maximum value is 1, this means that only 28% of the change in the dependent variable is attributable by the independent variables. Is this right? Shouldn't this be classified as poor/ fair?


title: statistics on proportion wise data

text: What are the usual statistical procedures used for data that is in proportions ? I have data about survival proportions of two species across 7 elevations .. What statistical tests could I run apart from regression on proportions across the 7 elevations ?


title: How to test homoscedasticity when the errors are DEPENDENT?

text: I have done a linear regression and plot the residuals. I noticed that the errors are **dependent** (autocorrelated), how could I test homoscedasticity of this series? I read that Breusch-Pagan test only work with independent errors. What test could I use to test if the variance is constant in this case? Thank you


title: Convergence in probability vs. almost sure convergence

text: I've never really grokked the difference between these two measures of convergence. (Or, in fact, any of the different types of convergence, but I mention these two in particular because of the Weak and Strong Laws of Large Numbers.) Sure, I can quote the definition of each and give an example where they differ, but I still don't quite get it. What's a good way to understand the difference? Why is the difference important? Is there a particularly memorable example where they differ?


title: What is the null hypothesis in the Mann-Whitney test?

text: Let $X_1$ be a random value from distribution 1 and let $X_2$ be a random value from distribution 2. I thought that the null hypothesis for the Mann- Whitney test was $P(X_1 < X_2) = P(X_2 < X_1)$. If I run simulations of the Mann-Whitney test on data from normal distributions with equal means and equal variances, with $\alpha=0.05$, I get Type I error rates which are very close to 0.05. However, if I make the variances unequal (but leave the means equal), the proportion of simulations in which the null hypothesis is rejected becomes larger than 0.05, which I didn't expect, since $P(X_1 < X_2) = P(X_2 < X_1)$ still holds. This happens when I use `wilcox.test` in R, regardless of whether I have `exact=TRUE`, `exact=FALSE, correct=TRUE`, or `exact=FALSE, correct=FALSE`. Is the null hypothesis something different from what I've written above, or is it just that the test is inaccurate in terms of Type I error if the variances are unequal?


title: Deciphering output from R2WinBUGS

text: I have been struggling to figure out how to identify which output the parameter came from. For e.g. in the following code I run the model for 2 iterations. At the end of it I obtain 2 outputs named `bugs.output[[1]], bugs.output[[2]]`. Now when I type `mean (p[,2])` it gives me the mean for `p[2]` from the second output. How do I specify that I need the mean for `p[2]` for the `1st` or the `ith` output? Here is how my output looks:                        mean       sd      2.5%       25%       50%       75%     97.5%     alpha      4.20338  1.14447   2.44975   3.39950   4.01550   4.82925   7.05405     p[1]       0.16691  0.04036   0.09974   0.13720   0.16445   0.19222   0.25308     p[2]       0.14286  0.03812   0.07665   0.11670   0.14020   0.16702   0.22630     p[3]       0.21238  0.04240   0.13740   0.18295   0.21050   0.23845   0.30360.......      And here is my code:               library("R2WinBUGS")     trial.data <- read.table("simuldatBB_6_30.csv", header=T)     p_true<- read.table("p_BB_6_30.csv",header=F)     bugs.output <- list()     for(i in 1:2){            nausea <- as.integer(trial.data[i,])            bugs.output[[i]] <- bugs(            data=list(nausea=nausea, N=63),            inits=list(                    list(theta=.300,mu=3),                    list(theta=.350, mu=3)                    ),            model.file="conj_nausea_script_2.txt",            parameters.to.save = c("alpha","p"),     n.chains=2, n.iter=12000, n.burnin=5000,     bugs.directory="E://AChaudhuri/winbugs14/WinBUGS14",     working.directory=NULL)     attach.bugs(bugs.output[[i]],overwrite=NA)}


title: Quality-price trade-off

text: Say there are 3 companies A, B and C. Each company has a quality rating from 0 to 100 and a price in USD.               Company    Quality    Price     A          80         7.9     B          70         8.0     C          75         8.1      How do I determine the best quality-price trade-off? What kind of analysis should I use?


title: Select best distance for feature selection

text: Suppose I have matrix $X \in R^{n \times m}$, where $n$ is the number of individuals and $m$ is the number of features and $X[i,j] \in \\{0,1\\}$; $1$ indicates that the individual $i$ has the feature $j$, and $0$ means not. Individuals only belong to one of the two classes. Given a distance metric $Dist$ (any distance), I want to select the best set of features that gives me the maximum $KLD$ between the two classes. Is there any **feature selection** algorithm that can **given a distance** select the best features? Ultimately, I would like to test a bunch of distances.


title: How to estimate loyalty to service?

text: Let's say there are two services $A$ and $B$ that solve same problem. We would like to know which one users more loyal to. Loyalty $L(X)$ is probability that user will use service $X$ when asked to solve problem. Well for every user we can estimate $L$ for service $A$ like this $L(A) = \frac{N(A)}{N(A) + N(B)}$ where $N(X)$ is number of uses of $X$. We can even smooth it with Laplacian smoothing and get $L'(A) = \frac{N(A) + 1}{N(A) + N(B) + 2}$ but it's not the case. In reality 20% of users have $N(A) + N(B) = 1$, then there are two types of $L(A)$: 0 and 1 both with variance like 0.5 =). Maybe there are known methods how to behave in such situations when there was almost none experiments conducted. We need to have better estimation on both $L$ and it's variance.


title: Interpreting multivariate tests in repeated measures

text: I have a within subjects design with 4 DVs and I ran a repeated measures analysis on SPSS. My DV is a physiological measurement and is not expected to be linear. My multivariate Wilks lambda is significant at p = .004 and my Mauchly test approaches significance at p = .054. However, my univariate tests of within subjects is non-significant (uncorrected p =.07 and greenhouse Geisser p = .082)   1. Given that my sphericity test is almost significant and given that I do not expect linearity in the data, is it reasonable to look at my multivariate test and ignore the non-significant univariate test?    2. Do I even have to look at univariate if the multivariate is significant? I.e., am I justified to look at my pairwise comparisons simply because the multivariate is significant?


title: High dimensional time series

text: I'm not sure what words I should look for. I have an under determined dataset of 8000 correlated variables (sales) over 12 months (ie 12 observations for each variable). And I basically want to predict the future. Where should I start? PCA? My question is, what are the techniques used to deal with lots of (correlated) variables and few observations in the case of time series. I'm looking for orientations (but a full solution is of course welcome!). I'm agnostic, so stats/econometrics/model-based are as fine as machine learning/AI/not-model-based, as long as they yield results that are useful and mean something. Repeatable/standardizable solutions are more than welcome.


title: Is this a problem for Survival analysis?

text: I have a dataset of individuals. Each individual has the same start time at which we begin observing them. There is also an end time for all individuals. Some individuals fail before they reach the end time and some individuals never fail and reach the end time (i.e. they succeed). This is a survival analysis problem in that I am trying to model a time to event (failure). Where my confusion comes in is with the individuals who succeed. I obviously cannot treat these individuals as observed failures, but I also cannot treat them as censored. This is because censorship implies that they were not observed to fail but they will fail at some point in the future (we just don't know when). This is not true in my case because if an individual does not fail up to the end time, it will never fail. So how do I deal with these individuals who never fail? Is this still a survival analysis problem? Here is an example. Say a group of athletes all begin training for the next olympic games. Some of them will end up getting injured from training and will not be able to compete. But if an athlete does not get injured until the date of the olympic games, then he/she would be able to compete (i.e. succeed). So these athletes that don't get injured and are able to compete are not censored, because whether they get injured after the olympics is irrelevant. All we care about is whether they get injured up to the date of the olympics


title: Logistic regression with longitudinal data

text: I have longitudinal data (one measurement being continous and the other ordinal) and I want to fit a model which takes into account the correlation between these two measures. I am using stata. Can anyone help me? Thanks in advance


title: Creating univariable smoothed scatterplot on logit scale using R

text: I am learning logistic regression modeling from the book _Applied Logistic Regression_ by Hosmer. I need to create a plot named "create univariable smoothed scatterplot on logit scale", something like this one (Figure 4.2 page 107): ![enter image description here](http://i.stack.imgur.com/1HRT2.gif) Can anyone help? Thanks. ### edit 01 Thanks for all the answers. I tried creating the plot, and discover that `smooth.spline` created the graph by using percentage (`nrow((df[df[,"dfree"]==1,]))/nrow(df)`), rather than logit (`log(nrow(df[df[,"dfree"]==1,])/nrow(df[df[,"dfree"]==0,]))`) when creating the y-axis. Yes, the graph looks similar, but I wonder if we can create an exact copy of that? Thanks. One more thing, when using logit, some of the logit value in the data.frame is `-Inf` and `Inf`, which is not allowed in `smooth.spline`. ### edit 02 I have further question concerning my edit: What if I manually dropped all `Inf` and `-Inf` in my data.frame and do `smooth.spline` on that data.frame? Is it appropriate?


title: Alternatives to svychisq() in survey

text: A time ago I posted this question. Now I have a 'little' continuation: I'm having some difficult with a few questions, getting this error message:               > svychisq(~ p13 + p14, amostra)     Error in solve.default(denom, numr) :        system is computationally singular: reciprocal condition number = 4.64501e-18      There is any alternative package/analisys? Or, at least, an explanation to my boss of why I can't analyse this data :P EDIT 2: Data and commands (sorry for the long `dput()`). p11 is the stratum, poulacao is the stratum population size, p13 and p14 are two questions.               require(survey)     dados_censo <- structure(list(p11 = c(16L, 11L, 12L, 12L, 12L, 12L, 12L, 12L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 8L, 8L, 8L, 8L, 8L, 8L, 12L,      12L, 12L, 12L, 12L, 12L, 12L, 12L, 26L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 9L, 9L, 9L, 9L,      9L, 9L, 9L, 9L, 17L, 17L, 17L, 17L, 22L, 22L, 22L, 22L, 22L,      22L, 22L, 17L, 17L, 17L, 18L, 11L, 12L, 18L, 18L, 18L, 18L, 16L,      6L, 6L, 6L, 6L, 6L, 17L, 18L, 17L, 27L, 27L, 27L, 6L, 6L, 9L,      18L, 18L, 11L, 2L, 2L, 11L, 11L, 11L, 11L, 27L, 27L, 27L, 27L,      16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 3L, 3L, 3L, 3L,      3L, 3L, 3L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 22L, 22L, 16L, 16L,      16L, 16L, 16L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 8L, 18L, 18L, 10L,      16L, 16L, 3L, 3L, 3L, 16L, 26L, 20L, 26L, 9L, 1L, 9L, 5L, 5L,      5L, 5L, 5L, 5L, 5L, 5L, 5L, 13L, 5L, 6L, 6L, 6L, 6L, 1L, 18L,      13L, 27L, 27L, 27L, 13L, 13L, 6L, 6L, 10L, 10L, 9L, 10L, 10L,      9L, 10L, 6L, 6L, 21L, 21L, 21L, 21L, 21L, 6L, 6L, 5L, 6L, 13L,      13L, 13L, 13L, 13L, 18L, 20L, 20L, 20L, 17L, 17L, 21L, 10L, 10L,      14L, 6L, 27L, 27L, 27L, 27L, 13L, 13L, 3L, 24L, 9L, 24L, 24L,      24L, 24L, 9L, 9L, 13L, 13L, 2L, 2L, 2L, 14L, 14L, 14L, 14L, 14L,      14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,      14L, 14L, 14L, 8L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,      11L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 9L, 13L, 19L, 19L, 19L, 8L, 8L, 8L,      19L, 13L, 13L, 13L, 13L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,      15L, 15L, 15L, 15L, 19L, 19L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 13L, 13L, 13L,      14L, 14L, 13L, 13L, 9L, 8L, 27L, 25L, 25L, 25L, 7L, 7L, 7L, 24L,      16L, 16L, 6L, 6L, 11L, 13L, 11L, 26L, 28L, 28L, 28L, 28L, 28L,      28L, 28L, 28L, 28L, 28L, 28L, 13L, 13L, 13L, 13L, 13L, 10L, 10L,      10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 13L, 13L,      13L, 14L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 2L, 14L, 3L, 27L, 27L,      27L, 27L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 4L, 19L, 19L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 19L, 19L, 19L,      19L, 19L, 26L, 19L, 19L, 19L, 26L, 26L, 26L, 26L, 26L, 6L, 26L,      26L, 27L, 7L, 7L, 19L, 27L, 18L, 18L, 18L, 18L, 18L, 26L, 26L,      20L, 26L, 19L, 19L, 19L, 19L, 25L, 18L, 18L, 18L, 18L, 18L, 18L,      18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,      18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 25L, 3L, 3L, 3L, 7L,      19L, 19L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 24L,      24L, 24L, 14L, 14L, 14L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 13L, 13L, 13L, 13L, 13L, 13L, 18L, 18L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 18L, 24L, 19L, 19L, 19L, 19L,      19L, 19L, 14L, 2L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L,      25L, 25L, 25L, 25L, 25L, 25L, 24L, 24L, 24L, 24L, 8L, 8L, 8L,      8L, 13L, 13L, 13L, 13L, 18L, 18L, 18L, 18L, 13L, 13L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 26L, 24L, 8L, 8L, 19L, 19L,      19L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 24L, 26L, 26L, 26L, 24L, 18L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 26L, 26L, 26L, 26L, 24L, 24L, 24L, 24L, 26L, 24L, 14L,      14L, 14L, 27L, 27L, 27L, 27L, 25L, 25L, 25L, 25L, 25L, 25L, 26L,      26L, 16L, 6L, 6L, 6L, 6L, 6L, 6L, 19L, 19L, 19L, 19L, 19L, 10L,      19L, 14L, 14L, 14L, 14L, 14L, 14L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 26L, 26L, 24L, 20L,      27L, 27L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 11L, 11L, 11L, 11L,      13L, 11L, 11L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 16L, 16L, 27L, 19L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 11L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 26L, 26L, 26L,      19L, 13L, 13L, 13L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 13L, 13L, 24L, 24L, 24L, 24L, 26L, 26L,      26L, 26L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L, 13L,      13L, 13L, 13L, 13L, 13L, 23L, 23L, 23L, 13L, 13L, 13L, 26L, 26L,      26L, 26L, 26L, 26L, 26L, 13L, 13L, 13L, 13L, 26L, 26L, 22L, 22L,      22L, 8L, 8L, 8L, 13L, 26L, 11L, 13L, 23L, 23L, 15L, 23L, 13L,      13L, 13L, 11L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 11L, 9L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      19L, 19L, 19L, 19L, 19L, 19L, 10L, 19L, 13L, 13L, 19L, 19L, 19L,      13L, 13L, 13L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 13L, 13L,      13L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 3L, 19L, 13L,      13L, 13L, 13L, 21L, 21L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 19L, 19L, 19L, 21L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 24L, 24L, 24L, 24L, 24L, 11L, 14L, 14L, 15L,      13L, 13L, 21L, 25L, 25L, 25L, 24L, 24L, 24L, 24L, 24L, 19L, 19L,      19L, 19L, 19L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 27L, 24L, 24L,      24L, 24L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 24L, 24L,      11L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 19L, 19L, 11L, 11L, 11L, 11L, 19L, 19L, 19L,      19L, 19L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 19L, 19L,      13L, 13L, 13L, 13L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 13L, 13L, 19L, 13L,      13L, 13L, 13L, 19L, 19L, 24L, 24L, 13L, 19L, 19L, 19L, 13L, 13L,      13L, 13L, 13L, 13L, 14L, 25L, 19L, 19L, 11L, 19L, 11L, 11L, 11L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 19L, 11L, 13L, 13L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 6L, 6L, 13L, 13L, 13L, 13L, 19L,      13L, 19L, 19L, 19L, 19L, 19L, 13L, 19L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 13L, 19L, 19L, 14L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 13L, 13L, 13L, 13L,      11L, 19L, 19L, 13L, 11L, 13L, 13L, 13L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 19L, 19L, 13L, 13L, 19L, 19L, 19L, 19L, 19L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 19L, 13L, 13L, 19L, 13L, 13L, 13L,      11L, 11L, 11L, 21L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 19L, 19L, 19L, 19L, 19L, 13L,      11L, 24L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 24L, 24L, 24L, 24L, 13L, 13L, 15L, 14L, 2L, 25L, 26L,      26L, 26L, 26L, 26L, 12L, 12L, 12L, 12L, 8L, 8L, 8L, 19L, 19L,      19L, 19L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 23L, 13L, 13L, 13L, 13L, 13L, 13L, 19L, 13L, 19L, 19L, 13L,      13L, 24L, 24L, 13L, 13L, 13L, 26L, 26L, 24L, 19L, 11L, 19L, 19L,      19L, 19L, 19L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,      13L, 13L, 19L, 19L, 19L, 19L, 19L, 24L, 13L, 24L, 24L, 19L, 13L,      11L, 11L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 13L, 19L, 19L, 13L,      13L, 13L, 13L, 13L, 9L, 13L, 13L, 13L, 13L, 13L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 13L, 13L, 13L, 19L, 19L, 13L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 18L, 18L, 18L, 18L, 18L, 19L, 19L,      26L, 26L, 26L, 26L, 26L, 25L, 25L, 25L, 13L, 13L, 13L, 11L, 11L,      11L, 11L, 9L, 9L, 9L, 2L, 13L, 11L, 19L, 18L, 19L, 13L, 13L,      13L, 19L, 19L, 18L, 25L, 13L, 13L, 13L, 13L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 13L, 13L, 13L, 13L, 13L, 19L, 19L, 19L, 19L, 19L,      8L, 8L, 8L, 8L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 11L, 11L, 11L, 19L, 25L, 25L, 25L, 25L, 25L, 25L,      25L, 25L, 25L, 11L, 19L, 11L, 11L, 18L, 26L, 26L, 26L, 26L, 26L,      26L, 26L, 26L, 26L, 11L, 13L, 13L, 13L, 18L, 13L, 13L, 13L, 13L,      13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 25L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 18L, 19L, 19L, 19L, 19L,      19L, 13L, 13L, 18L, 19L, 19L, 13L, 25L, 19L, 19L, 13L, 13L, 25L,      13L, 25L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 26L, 26L, 26L, 26L, 26L, 26L,      13L, 13L, 13L, 25L, 25L, 25L, 25L, 25L, 13L, 13L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L), p13 = c(4L,      4L, 5L, 4L, 5L, 7L, 2L, 5L, 5L, 7L, 7L, 7L, 6L, 4L, 6L, 6L, 7L,      4L, 1L, 1L, 7L, 8L, 8L, 7L, 6L, 6L, 2L, 7L, 6L, 5L, 1L, 1L, 4L,      5L, 8L, 5L, 7L, 6L, 6L, 7L, 7L, 7L, 8L, 3L, 1L, 6L, 6L, 6L, 4L,      7L, 6L, 6L, 8L, 4L, 8L, 8L, 8L, 8L, 6L, 7L, 6L, 6L, 8L, 7L, 7L,      8L, 6L, 8L, 7L, 6L, 8L, 8L, 6L, 7L, 6L, 7L, 5L, 6L, 8L, 4L, 5L,      7L, 6L, 5L, 1L, 4L, 7L, 2L, 8L, 2L, 1L, 2L, 2L, 2L, 1L, 7L, 4L,      1L, 2L, 6L, 7L, 8L, 6L, 1L, 3L, 5L, 6L, 6L, 1L, 5L, 6L, 5L, 1L,      7L, 6L, 5L, 5L, 5L, 7L, 5L, 5L, 6L, 7L, 7L, 3L, 4L, 1L, 4L, 5L,      7L, 8L, 4L, 3L, 8L, 1L, 5L, 1L, 8L, 1L, 1L, 8L, 2L, 6L, 5L, 2L,      1L, 1L, 4L, 2L, 1L, 8L, 8L, 4L, 6L, 7L, 7L, 1L, 8L, 4L, 7L, 5L,      7L, 6L, 1L, 5L, 5L, 7L, 7L, 6L, 4L, 8L, 6L, 8L, 6L, 1L, 5L, 7L,      6L, 4L, 4L, 6L, 3L, 7L, 8L, 7L, 6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 7L, 1L, 1L, 7L, 6L, 8L, 6L, 6L, 1L, 1L, 6L, 8L, 4L, 7L,      1L, 7L, 6L, 1L, 5L, 1L, 5L, 7L, 8L, 1L, 5L, 1L, 1L, 1L, 1L, 7L,      5L, 1L, 1L, 1L, 5L, 4L, 5L, 7L, 5L, 8L, 8L, 8L, 4L, 4L, 8L, 5L,      7L, 6L, 5L, 8L, 8L, 8L, 6L, 7L, 6L, 1L, 4L, 6L, 4L, 1L, 1L, 5L,      7L, 8L, 5L, 5L, 3L, 4L, 7L, 2L, 8L, 1L, 5L, 5L, 1L, 8L, 1L, 6L,      5L, 7L, 8L, 7L, 4L, 8L, 8L, 4L, 5L, 3L, 7L, 4L, 4L, 8L, 6L, 3L,      6L, 5L, 8L, 7L, 4L, 1L, 8L, 8L, 5L, 5L, 5L, 1L, 8L, 7L, 7L, 8L,      8L, 4L, 6L, 4L, 1L, 6L, 6L, 8L, 6L, 6L, 6L, 1L, 7L, 3L, 7L, 7L,      6L, 4L, 6L, 6L, 5L, 1L, 8L, 2L, 1L, 2L, 2L, 5L, 2L, 8L, 1L, 1L,      4L, 4L, 5L, 4L, 6L, 5L, 5L, 3L, 6L, 5L, 4L, 4L, 8L, 8L, 7L, 8L,      8L, 8L, 7L, 4L, 3L, 7L, 6L, 1L, 7L, 8L, 8L, 8L, 6L, 4L, 7L, 8L,      6L, 5L, 5L, 8L, 5L, 7L, 1L, 5L, 1L, 1L, 8L, 5L, 8L, 4L, 1L, 1L,      1L, 1L, 1L, 8L, 8L, 1L, 1L, 3L, 6L, 7L, 8L, 1L, 1L, 3L, 6L, 8L,      6L, 1L, 1L, 1L, 1L, 7L, 1L, 6L, 6L, 8L, 1L, 2L, 1L, 4L, 3L, 1L,      1L, 5L, 7L, 6L, 1L, 1L, 4L, 7L, 2L, 1L, 4L, 7L, 4L, 2L, 7L, 2L,      2L, 2L, 1L, 3L, 6L, 4L, 7L, 4L, 6L, 3L, 5L, 1L, 2L, 1L, 8L, 4L,      6L, 5L, 7L, 7L, 7L, 8L, 7L, 6L, 6L, 8L, 8L, 8L, 4L, 1L, 7L, 1L,      5L, 4L, 4L, 7L, 4L, 1L, 1L, 5L, 4L, 4L, 6L, 4L, 7L, 8L, 2L, 3L,      2L, 4L, 1L, 4L, 6L, 6L, 7L, 6L, 7L, 8L, 8L, 6L, 6L, 6L, 6L, 7L,      8L, 6L, 1L, 5L, 6L, 7L, 1L, 4L, 5L, 1L, 6L, 8L, 7L, 7L, 4L, 8L,      7L, 8L, 8L, 7L, 5L, 7L, 6L, 7L, 7L, 5L, 7L, 7L, 7L, 6L, 6L, 5L,      6L, 5L, 5L, 7L, 7L, 4L, 6L, 4L, 7L, 7L, 8L, 4L, 7L, 7L, 4L, 5L,      3L, 2L, 1L, 1L, 6L, 5L, 2L, 1L, 5L, 7L, 3L, 3L, 5L, 6L, 5L, 8L,      5L, 4L, 1L, 5L, 6L, 6L, 7L, 1L, 7L, 6L, 7L, 6L, 5L, 7L, 7L, 5L,      5L, 2L, 5L, 1L, 1L, 7L, 7L, 7L, 6L, 8L, 7L, 5L, 5L, 7L, 8L, 6L,      6L, 7L, 8L, 4L, 4L, 4L, 2L, 4L, 5L, 4L, 5L, 8L, 8L, 5L, 5L, 8L,      8L, 5L, 8L, 8L, 3L, 4L, 4L, 4L, 5L, 1L, 4L, 7L, 7L, 6L, 8L, 5L,      8L, 5L, 6L, 4L, 8L, 5L, 1L, 4L, 8L, 8L, 1L, 4L, 5L, 5L, 5L, 5L,      6L, 8L, 5L, 5L, 8L, 6L, 7L, 5L, 7L, 1L, 3L, 4L, 4L, 6L, 5L, 2L,      5L, 8L, 4L, 1L, 5L, 5L, 5L, 4L, 1L, 8L, 6L, 6L, 8L, 6L, 8L, 7L,      2L, 6L, 8L, 8L, 8L, 7L, 7L, 6L, 8L, 6L, 6L, 8L, 2L, 6L, 8L, 5L,      7L, 6L, 5L, 6L, 6L, 5L, 4L, 7L, 5L, 4L, 4L, 5L, 6L, 4L, 7L, 7L,      6L, 6L, 5L, 6L, 2L, 2L, 1L, 4L, 1L, 1L, 4L, 4L, 6L, 1L, 1L, 4L,      5L, 7L, 4L, 4L, 5L, 6L, 2L, 8L, 3L, 1L, 5L, 1L, 1L, 8L, 6L, 2L,      7L, 6L, 8L, 7L, 6L, 3L, 8L, 4L, 7L, 4L, 1L, 7L, 6L, 4L, 3L, 3L,      6L, 4L, 1L, 4L, 1L, 1L, 4L, 4L, 5L, 6L, 6L, 4L, 6L, 1L, 1L, 5L,      5L, 7L, 6L, 1L, 2L, 4L, 3L, 4L, 1L, 4L, 6L, 6L, 7L, 7L, 4L, 7L,      8L, 2L, 4L, 1L, 3L, 4L, 7L, 4L, 5L, 1L, 6L, 7L, 1L, 3L, 1L, 2L,      4L, 1L, 6L, 4L, 2L, 3L, 2L, 4L, 3L, 1L, 2L, 1L, 5L, 6L, 2L, 8L,      6L, 1L, 6L, 6L, 5L, 2L, 6L, 4L, 8L, 4L, 7L, 4L, 8L, 8L, 6L, 6L,      7L, 7L, 7L, 6L, 3L, 4L, 2L, 5L, 8L, 8L, 6L, 5L, 5L, 5L, 8L, 5L,      6L, 3L, 6L, 2L, 1L, 4L, 8L, 7L, 7L, 7L, 7L, 8L, 5L, 8L, 6L, 6L,      8L, 4L, 6L, 7L, 6L, 5L, 1L, 6L, 8L, 8L, 7L, 4L, 8L, 6L, 7L, 7L,      1L, 7L, 8L, 5L, 3L, 6L, 3L, 6L, 5L, 4L, 1L, 7L, 1L, 5L, 2L, 6L,      4L, 6L, 7L, 8L, 1L, 1L, 8L, 6L, 6L, 1L, 5L, 3L, 6L, 6L, 7L, 2L,      8L, 6L, 6L, 6L, 1L, 8L, 6L, 4L, 8L, 1L, 1L, 5L, 2L, 3L, 7L, 7L,      1L, 4L, 6L, 1L, 4L, 1L, 5L, 6L, 5L, 4L, 3L, 4L, 6L, 1L, 3L, 2L,      6L, 4L, 6L, 4L, 5L, 7L, 1L, 2L, 5L, 6L, 6L, 6L, 2L, 5L, 4L, 3L,      2L, 6L, 2L, 1L, 6L, 7L, 6L, 5L, 5L, 5L, 4L, 1L, 1L, 4L, 4L, 1L,      6L, 5L, 5L, 4L, 1L, 1L, 1L, 1L, 1L, 6L, 4L, 7L, 4L, 2L, 6L, 1L,      1L, 1L, 8L, 4L, 1L, 4L, 5L, 1L, 4L, 5L, 1L, 2L, 5L, 6L, 1L, 5L,      7L, 8L, 4L, 2L, 1L, 2L, 1L, 3L, 8L, 7L, 8L, 7L, 4L, 4L, 1L, 4L,      8L, 1L, 2L, 7L, 3L, 5L, 5L, 6L, 6L, 5L, 4L, 8L, 5L, 3L, 5L, 4L,      4L, 6L, 5L, 7L, 1L, 8L, 4L, 7L, 6L, 8L, 1L, 7L, 5L, 4L, 7L, 5L,      6L, 5L, 6L, 4L, 4L, 5L, 7L, 6L, 3L, 6L, 5L, 7L, 4L, 7L, 4L, 1L,      1L, 8L, 6L, 6L, 4L, 8L, 7L, 7L, 2L, 3L, 1L, 6L, 4L, 7L, 3L, 4L,      1L, 1L, 6L, 5L, 4L, 7L, 1L, 5L, 5L, 1L, 1L, 4L, 6L, 1L, 5L, 6L,      3L, 1L, 4L, 8L, 6L, 4L, 5L, 3L, 6L, 5L, 4L, 3L, 6L, 5L, 1L, 5L,      6L, 8L, 4L, 6L, 2L, 5L, 6L, 4L, 1L, 5L, 1L, 1L, 7L, 3L, 6L, 2L,      1L, 5L, 6L, 1L, 7L, 4L, 5L, 7L, 4L, 6L, 2L, 2L, 7L, 1L, 7L, 2L,      2L, 8L, 7L, 4L, 4L, 6L, 1L, 2L, 6L, 2L, 8L, 8L, 8L, 7L, 4L, 5L,      2L, 4L, 3L, 1L, 1L, 2L, 6L, 6L, 5L, 3L, 8L, 3L, 6L, 1L, 3L, 5L,      8L, 4L, 4L, 7L, 8L, 5L, 4L, 5L, 7L, 5L, 3L, 6L, 1L, 3L, 4L, 3L,      6L, 4L, 6L, 8L, 6L, 6L, 7L, 8L, 7L, 5L, 5L, 6L, 6L, 8L, 7L, 1L,      5L, 6L, 6L, 4L, 6L, 6L, 1L, 6L, 1L, 6L, 5L, 7L, 3L, 5L, 1L, 4L,      5L, 7L, 3L, 6L, 1L, 4L, 7L, 1L, 5L, 5L, 5L, 7L, 7L, 4L, 1L, 3L,      4L, 6L, 7L, 2L, 6L, 7L, 6L, 5L, 4L, 2L, 3L, 3L, 2L, 5L, 1L, 2L,      3L, 8L, 5L, 7L, 4L, 6L, 5L, 8L, 4L, 4L, 6L, 1L, 6L, 6L, 4L, 7L,      7L, 2L, 7L, 5L, 6L, 8L, 8L, 8L, 6L, 7L, 5L, 7L, 2L, 5L, 4L, 4L,      4L, 1L, 1L, 8L, 6L, 1L, 4L, 7L, 6L, 5L, 4L, 1L, 6L, 6L, 6L, 6L,      6L, 1L, 1L, 1L, 4L, 7L, 4L, 5L, 4L, 1L, 3L, 6L, 7L, 7L, 3L, 4L,      8L, 6L, 5L, 4L, 8L, 6L, 5L, 6L, 6L, 4L, 7L, 6L, 1L, 1L, 4L, 2L,      7L, 7L, 3L, 5L, 5L, 8L, 7L, 2L, 5L, 2L, 5L, 4L, 5L, 4L, 8L, 1L,      6L, 4L, 1L, 1L, 1L, 5L, 3L, 5L, 1L, 1L, 1L, 4L, 4L, 6L, 6L, 4L,      1L, 3L, 1L, 5L, 7L, 1L, 6L, 5L, 1L, 5L, 5L, 2L, 4L, 1L, 7L, 3L,      2L, 2L, 6L, 6L, 8L, 5L, 4L, 1L, 1L, 6L, 1L, 4L, 6L, 6L, 6L, 8L,      6L, 6L, 5L, 1L, 3L, 1L, 3L, 5L, 6L, 6L, 5L, 6L, 2L, 6L, 7L, 1L,      6L, 7L, 6L, 1L, 8L, 4L, 8L, 7L, 6L, 7L, 6L, 8L, 7L, 6L, 8L, 7L,      6L, 5L, 6L, 3L, 4L, 4L, 1L, 6L, 7L, 5L, 5L, 7L, 7L, 5L, 3L, 8L,      6L, 5L, 8L, 5L, 6L, 6L, 6L, 2L, 6L, 6L, 4L, 5L, 1L, 7L, 6L, 6L,      7L, 6L, 8L, 2L, 1L, 6L, 6L, 7L, 4L, 6L, 1L, 4L, 6L, 7L, 5L, 4L,      1L, 4L, 5L, 4L, 1L, 6L, 6L, 6L, 4L, 6L, 1L, 6L, 4L, 6L, 1L, 4L,      8L, 6L, 7L, 6L, 1L, 2L, 4L, 1L, 1L, 4L, 6L, 6L, 1L, 3L, 4L, 1L,      7L, 2L, 6L, 2L, 5L, 4L, 7L, 6L, 4L, 4L, 5L, 5L, 6L, 4L, 4L, 7L,      4L, 5L, 7L, 6L, 6L, 6L, 5L, 6L, 4L, 4L, 1L, 2L, 4L, 6L, 4L, 8L,      6L, 5L, 4L, 4L, 6L, 4L, 5L, 4L, 4L, 5L, 4L, 4L, 6L, 4L, 5L, 4L,      7L, 7L, 6L, 6L, 6L, 2L, 6L, 6L, 6L, 7L, 7L, 6L, 4L, 7L, 4L, 6L,      5L, 6L, 7L, 4L, 6L, 8L, 8L, 4L, 4L, 6L, 3L, 7L, 4L, 7L, 5L, 7L,      7L, 1L, 1L, 6L, 6L, 7L, 4L, 4L, 5L, 1L, 1L, 1L, 6L, 1L, 3L, 4L,      6L, 6L, 7L, 4L, 5L, 6L, 4L, 2L, 1L, 4L, 1L, 7L, 6L, 1L, 4L, 2L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 6L, 4L, 6L, 7L, 7L, 8L, 4L, 6L,      7L, 8L, 5L, 8L, 3L, 5L, 7L, 4L, 7L, 4L, 1L, 5L, 7L, 8L, 5L, 8L,      6L, 6L, 7L, 6L, 7L, 8L, 5L, 7L, 4L, 7L, 7L, 5L, 8L, 8L, 5L, 1L,      6L, 6L, 3L, 5L, 1L, 5L, 6L, 4L, 8L, 1L, 7L, 7L, 4L, 8L, 1L, 5L,      7L, 4L, 2L, 5L, 1L, 4L, 8L, 6L, 7L, 5L, 7L, 7L, 6L, 5L, 1L, 5L,      5L, 7L, 6L, 1L, 5L, 8L, 5L, 5L, 5L, 7L, 5L, 1L, 5L, 3L, 6L, 4L,      5L, 6L, 8L, 8L, 5L, 5L, 6L, 1L, 4L, 4L, 6L, 7L, 7L, 7L, 6L, 8L,      5L, 7L, 6L, 8L, 4L, 5L, 5L, 7L, 6L, 2L, 5L, 1L, 3L, 5L, 6L, 8L,      3L, 3L, 4L, 3L, 6L, 3L, 4L), p14 = c(7L, 9L, 5L, 7L, 2L, 2L,      6L, 7L, 6L, 2L, 3L, 2L, 2L, 4L, 2L, 3L, 2L, 7L, 7L, 8L, 2L, 7L,      5L, 4L, 3L, 5L, 5L, 6L, 5L, 8L, 7L, 6L, 7L, 2L, 2L, 3L, 7L, 5L,      3L, 1L, 8L, 2L, 7L, 6L, 6L, 7L, 7L, 7L, 8L, 2L, 4L, 5L, 6L, 7L,      7L, 5L, 2L, 5L, 3L, 2L, 8L, 5L, 5L, 2L, 2L, 2L, 4L, 3L, 3L, 10L,      2L, 5L, 3L, 3L, 5L, 2L, 7L, 8L, 2L, 7L, 7L, 7L, 8L, 7L, 8L, 7L,      5L, 6L, 5L, 6L, 5L, 7L, 9L, 7L, 5L, 7L, 8L, 6L, 8L, 3L, 7L, 6L,      7L, 6L, 3L, 5L, 6L, 3L, 7L, 8L, 6L, 6L, 7L, 3L, 3L, 7L, 8L, 5L,      5L, 6L, 2L, 5L, 5L, 5L, 7L, 5L, 5L, 9L, 8L, 4L, 7L, 6L, 8L, 8L,      6L, 2L, 7L, 7L, 6L, 6L, 5L, 5L, 6L, 8L, 7L, 9L, 5L, 8L, 3L, 6L,      3L, 2L, 9L, 2L, 2L, 2L, 7L, 2L, 7L, 8L, 5L, 6L, 8L, 7L, 7L, 7L,      5L, 7L, 8L, 8L, 2L, 5L, 3L, 5L, 6L, 4L, 2L, 7L, 7L, 8L, 5L, 9L,      5L, 5L, 6L, 5L, 6L, 5L, 7L, 6L, 7L, 7L, 6L, 8L, 6L, 5L, 9L, 9L,      5L, 6L, 2L, 2L, 7L, 6L, 9L, 3L, 2L, 8L, 8L, 7L, 4L, 6L, 6L, 7L,      5L, 5L, 7L, 7L, 6L, 7L, 6L, 7L, 10L, 5L, 7L, 9L, 6L, 6L, 7L,      5L, 7L, 7L, 5L, 3L, 9L, 5L, 5L, 7L, 7L, 6L, 7L, 2L, 7L, 7L, 3L,      2L, 7L, 5L, 3L, 7L, 8L, 8L, 4L, 8L, 6L, 7L, 6L, 2L, 8L, 7L, 9L,      9L, 7L, 5L, 6L, 6L, 7L, 7L, 7L, 5L, 7L, 7L, 7L, 7L, 2L, 5L, 7L,      7L, 5L, 6L, 8L, 7L, 8L, 2L, 7L, 8L, 5L, 5L, 5L, 2L, 5L, 7L, 7L,      7L, 6L, 5L, 5L, 5L, 6L, 5L, 7L, 5L, 7L, 2L, 3L, 3L, 6L, 7L, 7L,      6L, 3L, 3L, 5L, 3L, 7L, 8L, 9L, 5L, 5L, 2L, 4L, 3L, 8L, 8L, 2L,      5L, 7L, 5L, 7L, 8L, 8L, 8L, 7L, 6L, 7L, 6L, 7L, 6L, 7L, 9L, 9L,      4L, 3L, 3L, 8L, 2L, 7L, 8L, 5L, 3L, 7L, 5L, 9L, 7L, 5L, 1L, 3L,      5L, 8L, 6L, 4L, 6L, 4L, 4L, 6L, 5L, 7L, 2L, 7L, 6L, 5L, 7L, 8L,      8L, 8L, 8L, 8L, 8L, 7L, 2L, 5L, 10L, 7L, 6L, 7L, 7L, 5L, 8L,      4L, 5L, 6L, 6L, 4L, 8L, 5L, 5L, 7L, 6L, 7L, 7L, 7L, 4L, 5L, 8L,      7L, 7L, 6L, 6L, 8L, 4L, 1L, 7L, 5L, 7L, 10L, 6L, 5L, 7L, 8L,      5L, 5L, 6L, 6L, 5L, 3L, 7L, 9L, 9L, 7L, 8L, 7L, 8L, 7L, 7L, 8L,      7L, 6L, 6L, 6L, 6L, 7L, 9L, 7L, 8L, 6L, 7L, 7L, 3L, 3L, 4L, 2L,      2L, 2L, 3L, 2L, 2L, 7L, 1L, 2L, 3L, 2L, 7L, 5L, 5L, 7L, 7L, 8L,      7L, 3L, 8L, 5L, 7L, 7L, 8L, 5L, 5L, 9L, 5L, 7L, 9L, 7L, 7L, 7L,      2L, 7L, 2L, 5L, 3L, 6L, 2L, 8L, 5L, 2L, 4L, 3L, 5L, 2L, 5L, 6L,      9L, 7L, 5L, 7L, 7L, 7L, 7L, 7L, 2L, 2L, 2L, 2L, 8L, 8L, 8L, 8L,      7L, 2L, 7L, 5L, 3L, 2L, 2L, 4L, 3L, 1L, 5L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 4L, 4L, 4L, 3L, 3L, 4L, 4L, 3L, 4L, 4L, 7L, 7L, 8L, 7L,      7L, 6L, 3L, 9L, 8L, 8L, 8L, 8L, 10L, 9L, 5L, 7L, 7L, 5L, 5L,      5L, 7L, 7L, 7L, 7L, 5L, 7L, 4L, 5L, 4L, 4L, 7L, 3L, 4L, 7L, 8L,      8L, 7L, 7L, 7L, 3L, 4L, 7L, 5L, 1L, 6L, 5L, 3L, 6L, 3L, 5L, 7L,      6L, 5L, 5L, 5L, 8L, 7L, 8L, 8L, 7L, 8L, 3L, 3L, 2L, 5L, 8L, 8L,      8L, 5L, 7L, 8L, 8L, 8L, 8L, 8L, 7L, 7L, 4L, 3L, 4L, 2L, 3L, 7L,      2L, 8L, 7L, 8L, 8L, 8L, 7L, 4L, 5L, 8L, 8L, 10L, 5L, 5L, 5L,      7L, 6L, 7L, 8L, 5L, 9L, 3L, 4L, 2L, 8L, 8L, 5L, 8L, 8L, 8L, 8L,      8L, 7L, 8L, 8L, 7L, 8L, 8L, 8L, 9L, 3L, 8L, 3L, 3L, 2L, 3L, 4L,      6L, 5L, 3L, 5L, 3L, 5L, 5L, 2L, 7L, 5L, 2L, 5L, 7L, 1L, 1L, 3L,      7L, 7L, 6L, 7L, 7L, 3L, 8L, 4L, 7L, 7L, 6L, 8L, 8L, 8L, 3L, 2L,      6L, 5L, 5L, 4L, 8L, 7L, 6L, 8L, 6L, 8L, 7L, 8L, 9L, 5L, 6L, 7L,      7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 7L, 8L, 7L, 5L, 3L, 5L, 6L,      1L, 6L, 1L, 2L, 2L, 5L, 5L, 8L, 8L, 6L, 6L, 5L, 8L, 5L, 6L, 6L,      8L, 9L, 5L, 7L, 6L, 6L, 7L, 10L, 8L, 8L, 8L, 5L, 8L, 8L, 8L,      5L, 8L, 2L, 8L, 5L, 8L, 8L, 8L, 5L, 8L, 6L, 4L, 5L, 5L, 3L, 8L,      3L, 3L, 8L, 7L, 8L, 7L, 7L, 8L, 8L, 7L, 8L, 7L, 8L, 8L, 8L, 8L,      8L, 8L, 6L, 6L, 5L, 8L, 8L, 9L, 9L, 8L, 8L, 9L, 7L, 8L, 8L, 7L,      7L, 9L, 7L, 8L, 4L, 5L, 6L, 7L, 5L, 6L, 8L, 7L, 8L, 1L, 5L, 6L,      4L, 3L, 5L, 2L, 7L, 8L, 8L, 8L, 7L, 7L, 7L, 4L, 5L, 5L, 4L, 5L,      8L, 7L, 7L, 7L, 6L, 6L, 8L, 4L, 4L, 4L, 3L, 7L, 6L, 7L, 4L, 7L,      5L, 1L, 8L, 8L, 7L, 5L, 4L, 6L, 1L, 4L, 2L, 2L, 7L, 7L, 5L, 2L,      2L, 7L, 7L, 2L, 7L, 6L, 8L, 8L, 7L, 8L, 8L, 8L, 5L, 8L, 7L, 8L,      8L, 8L, 3L, 2L, 5L, 2L, 3L, 2L, 3L, 7L, 8L, 8L, 6L, 2L, 2L, 5L,      9L, 5L, 2L, 5L, 5L, 6L, 3L, 6L, 8L, 7L, 7L, 8L, 7L, 10L, 9L,      7L, 2L, 5L, 8L, 9L, 5L, 8L, 5L, 6L, 5L, 8L, 8L, 5L, 5L, 6L, 7L,      8L, 8L, 9L, 7L, 7L, 9L, 5L, 5L, 7L, 3L, 8L, 8L, 8L, 10L, 9L,      5L, 5L, 8L, 9L, 8L, 8L, 7L, 5L, 7L, 7L, 6L, 7L, 7L, 8L, 9L, 9L,      8L, 7L, 7L, 5L, 7L, 7L, 9L, 6L, 8L, 8L, 8L, 8L, 10L, 7L, 7L,      8L, 9L, 7L, 8L, 6L, 7L, 7L, 5L, 5L, 8L, 8L, 8L, 7L, 9L, 8L, 8L,      8L, 8L, 7L, 6L, 8L, 9L, 8L, 8L, 8L, 7L, 7L, 8L, 7L, 8L, 7L, 8L,      8L, 8L, 8L, 9L, 7L, 5L, 9L, 7L, 10L, 2L, 6L, 8L, 8L, 8L, 5L,      1L, 4L, 6L, 5L, 7L, 7L, 7L, 8L, 8L, 5L, 7L, 7L, 8L, 5L, 1L, 5L,      3L, 6L, 5L, 3L, 9L, 10L, 8L, 6L, 8L, 8L, 5L, 8L, 3L, 8L, 3L,      5L, 3L, 5L, 3L, 9L, 6L, 8L, 5L, 7L, 5L, 7L, 8L, 3L, 3L, 7L, 7L,      8L, 8L, 8L, 7L, 6L, 6L, 7L, 6L, 7L, 4L, 8L, 10L, 7L, 5L, 8L,      9L, 8L, 7L, 4L, 5L, 6L, 5L, 8L, 6L, 5L, 8L, 5L, 5L, 8L, 7L, 2L,      9L, 8L, 7L, 7L, 7L, 8L, 8L, 9L, 8L, 8L, 8L, 9L, 8L, 3L, 8L, 6L,      4L, 5L, 5L, 8L, 8L, 8L, 8L, 8L, 6L, 4L, 8L, 2L, 5L, 5L, 8L, 8L,      4L, 8L, 8L, 5L, 5L, 10L, 4L, 5L, 7L, 3L, 8L, 9L, 3L, 7L, 7L,      8L, 8L, 7L, 7L, 3L, 6L, 8L, 8L, 8L, 8L, 8L, 7L, 9L, 8L, 8L, 8L,      8L, 7L, 4L, 7L, 5L, 8L, 5L, 6L, 2L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,      8L, 8L, 9L, 8L, 8L, 8L, 8L, 4L, 4L, 8L, 8L, 6L, 5L, 5L, 5L, 8L,      9L, 9L, 8L, 7L, 7L, 1L, 4L, 5L, 3L, 5L, 1L, 5L, 8L, 8L, 8L, 8L,      9L, 8L, 8L, 10L, 8L, 5L, 7L, 6L, 7L, 9L, 10L, 9L, 7L, 8L, 9L,      9L, 8L, 7L, 6L, 2L, 8L, 8L, 9L, 4L, 9L, 8L, 8L, 8L, 10L, 8L,      7L, 7L, 9L, 8L, 8L, 8L, 7L, 7L, 8L, 8L, 7L, 8L, 9L, 8L, 8L, 8L,      8L, 8L, 8L, 8L, 9L, 8L, 8L, 7L, 5L, 8L, 3L, 8L, 3L, 3L, 3L, 3L,      4L, 4L, 5L, 2L, 5L, 2L, 6L, 7L, 5L, 7L, 5L, 6L, 6L, 6L, 7L, 6L,      5L, 4L, 7L, 7L, 5L, 5L, 8L, 5L, 4L, 9L, 9L, 8L, 8L, 8L, 5L, 3L,      8L, 7L, 8L, 8L, 9L, 3L, 2L, 4L, 2L, 4L, 3L, 5L, 7L, 8L, 7L, 7L,      4L, 5L, 8L, 5L, 2L, 8L, 8L, 8L, 8L, 9L, 3L, 8L, 8L, 7L, 7L, 8L,      10L, 8L, 5L, 6L, 3L, 7L, 5L, 8L, 8L, 7L, 5L, 8L, 8L, 9L, 8L,      8L, 7L, 5L, 3L, 6L, 6L, 8L, 8L, 8L, 8L, 6L, 7L, 8L, 7L, 7L, 2L,      7L, 8L, 7L, 7L, 5L, 6L, 7L, 8L, 8L, 8L, 8L, 7L, 8L, 10L, 9L,      8L, 8L, 7L, 8L, 8L, 5L, 8L, 5L, 8L, 8L, 2L, 7L, 3L, 3L, 5L, 8L,      8L, 9L, 8L, 7L, 8L, 8L, 8L, 9L, 5L, 7L, 7L, 6L, 4L, 9L, 4L, 7L,      4L, 6L, 6L, 2L, 2L, 2L, 2L, 7L, 4L, 7L, 7L, 8L, 2L, 2L, 2L, 9L,      8L, 8L, 7L, 5L, 3L, 6L, 5L, 5L, 2L, 5L, 5L, 3L, 2L, 3L, 3L, 3L,      3L, 8L, 8L, 8L, 5L, 5L, 7L, 8L, 6L, 6L, 4L, 8L, 8L, 9L, 7L, 7L,      9L, 5L, 5L, 7L, 7L, 5L, 8L, 10L, 10L, 6L, 5L, 9L, 7L, 9L, 8L,      8L, 6L, 8L, 7L, 2L, 9L, 5L, 7L, 8L, 6L, 6L, 7L, 9L, 8L, 8L, 5L,      6L, 5L, 9L, 8L, 8L, 5L, 8L, 9L, 8L, 8L, 8L, 5L, 8L, 7L, 6L, 4L,      8L, 7L, 6L, 3L, 3L, 9L, 8L, 5L, 8L, 7L, 9L, 8L, 8L, 7L, 5L, 10L,      6L, 8L, 8L, 7L, 5L, 7L, 5L, 5L, 7L, 8L, 8L, 7L, 8L, 8L, 9L, 5L,      7L, 7L, 7L, 5L, 5L, 5L, 5L, 5L, 5L, 7L, 5L, 5L, 6L, 8L, 5L, 5L,      5L, 3L, 10L, 8L, 5L, 2L, 3L, 2L, 2L, 7L, 7L, 7L, 2L, 9L, 8L,      6L, 5L, 5L, 5L, 8L, 5L, 8L, 6L, 8L, 5L, 8L, 2L, 8L, 2L, 8L, 9L,      10L, 7L, 5L, 2L, 7L, 7L, 8L, 8L, 5L, 6L, 8L, 5L, 5L, 7L, 8L,      8L, 8L, 8L, 8L, 6L, 8L, 8L, 8L, 10L, 9L, 3L, 2L, 8L, 5L, 8L,      3L, 8L, 8L, 7L, 5L, 5L, 7L, 7L, 8L, 8L, 7L, 5L, 3L, 4L, 8L, 7L,      7L, 7L, 7L, 7L, 7L, 5L, 6L, 7L, 3L, 8L, 5L, 6L, 5L, 5L, 3L, 3L,      4L, 2L, 2L, 5L, 3L, 2L, 5L, 8L, 8L, 8L, 4L, 8L, 5L, 9L, 2L, 7L,      8L, 8L, 8L, 8L, 5L, 8L, 8L, 5L, 8L, 8L, 8L, 5L, 8L, 7L, 8L, 7L,      5L, 8L, 8L, 8L, 7L, 8L, 7L, 8L, 5L, 5L, 8L, 8L, 5L, 7L, 8L, 5L,      4L, 8L, 8L, 6L, 6L, 4L, 6L, 7L, 7L, 8L, 8L, 8L, 8L, 10L, 10L,      8L, 3L, 3L, 8L, 5L, 5L, 5L, 8L, 9L, 5L, 5L, 4L, 2L, 4L, 2L, 3L,      2L, 5L, 7L, 7L, 7L, 6L, 6L, 7L, 7L, 7L, 7L, 6L, 5L, 8L, 5L, 8L,      7L, 10L, 5L, 7L, 7L, 6L, 8L, 8L)), .Names = c("p11", "p13", "p14"     ), class = "data.frame", row.names = c(NA, -1768L))          populacao <- structure(list(p11 = 1:28, Populacao = c(225L, 115L, 119L, 9L,      36L, 310L, 8L, 568L, 230L, 193L, 569L, 145L, 4781L, 444L, 40L,      193L, 11L, 672L, 4417L, 101L, 21L, 15L, 140L, 2061L, 466L, 1204L,      302L, 87L)), .Names = c("p11", "Populacao"), row.names = c(NA,      -28L), class = "data.frame")          dados_censo <- merge(dados_censo,populacao)     amostra <- svydesign(id=~1,strata=~p11,  data=dados_censo, fpc=~Populacao)     svytable(~ p13 + p14, amostra)     svychisq(~ p13 + p14, amostra)


title: Gaussian process scale targets

text: I am currently playing around with Gaussian process regression. I discovered some confusing facts that I would like to be answered:   1. When I optimize the parameters like the length and error terms for the covariance function I get different results for scaled targets like $y \times 100$. This leads to the next thing:   2. When I use very small targets between $0$ and $0.1$ the whole process seems to perform very poor. Can somebody explain these behaviors to me or maybe can send a link where those problems are discussed?


title: Identification of confounder in a logistic regression model example in "Applied Logistic Regression"

text: I am reading Hosmer's _Applied Logistic Regression_ , and I am a bit stuck in chapter 3, when taking about interaction and confounders. In page 77, it states the following: > Using the estimated coefficient for LWD in model 1 we estimate the odds > ratio as $\exp(1.054)=2.87$. The results shown in Table 3.14 indicate that > AGE is not a strong confounder, $\Delta(\hat{\beta})\%=4.2$, but it does > interact with LWD, $p=0.076$. I am OK with the $\exp(1.054)=2.87$, but I can't figure out how he calculate the $\Delta(\hat{\beta})\%=4.2$. Can anyone help? Thanks.


title: Distance correlation versus mutual information

text: I've worked with the mutual information for some time. But I found a very recent measure in the "correlation world" that can also be used to measure distribution independence, the so called "distance correlation" ( also termed Brownian correlation): http://en.wikipedia.org/wiki/Brownian_covariance#Definitions . I checked the papers where this measure is introduced, but without finding any allusion to the mutual information. So, my questions are:   * Do they solve exactly the same problem? If not, how the problems are different?   * And if the previous question can be answered on the positive, what are the advantages of using one or the other? Thanks beforehand.


title: Including the interaction but not the main effects in a model

text: Is it ever valid to include a two-way interaction in a model without including the main effects? What if your hypothesis is only about the interaction, do you still need to include the main effects?


title: To exclude a main effect that is not of interest in ANCOVA, or not?

text: > **Possible Duplicate:**   >  Including the interaction but not the main effects in a model I have an experimental design with pretest and posttest. Two groups, experimental and control group. I am interested in finding whether there will be a difference between experimental and control group at the posttest, controlling for pretest. Therefore, I run ANCOVA with pretest as covariate, posttest as DV, and group as IV. I am also interested to see if there will be a gender differences for the experimental effect. Thus, I am interested in looking at the interaction term group*gender. I am not interested in looking at gender main effect alone, as it doesn't mean much to this experimental design. I customized the model and excluded the main effect of gender, and the results showed only the main effect of group and the interaction effect of group*gender. Is this okay, or do I need to run a full model including gender main effect? Thank you.


title: How to do regression with known correlations among the errors?

text: I was wondering if there is anything you can do when you have a regression problem: $$\begin{cases} Y_t = \beta_1x_t + \beta_0 + \varepsilon_t \\\ \left(\varepsilon_1,\ldots,\varepsilon_n\right)\sim\mathcal{N}_n\left(\mathbf{0},D \right) \\\ 1\leq t \leq n \end{cases}$$ and you know $D$, the covariance matrix of the errors. I know weighted least squares when $D$ is diagonal. But what happens if it isn't? Is there any well established estimator for $(\beta_0,\beta_1)$ ?


title: multi variable probability

text: I have a fair grip of probability and the normal distribution, but am looking at a problem associated with multiple variables, specifically magnitude and interval. I have a list of sizes and the time interval between them. The list I have been given is considered "ok". What I am looking for is the series of terms I should look up in order to that I can begin research into the solution to the problem of being given a new value and interval since the last, what is the probability that that value is also "ok". The second stage which I have no idea how to start is where there is a lack of event. I guess a field that would be closest to what I'm playing with is earthquakes where there are magnitude and interval (aftershocks notwithstanding) but my googling skills have not led me to what I need. Your thoughts would be most appreciated.


title: What method can be used to detect seasonality in data?

text: I want to detect seasonality in data that I receive. There are some methods that I have found like the seasonal subseries plot and the autocorrelation plot but the thing is I don't understand how to read the graph, could anyone help? The other thing is, are there other methods to detect seasonality with or without the final result in graph?


title: Test to compare user interfaces

text: For my thesis I will need to do an experimental user study to compare two user interfaces which both can perform the same tasks. They differ in how they support you in achieving your task. I'm not that familiar with statistics so please correct me if any of the following assumptions I already made are incorrect. I'm pretty sure that I'll need to use a **balanced within-subjects design**. Reading through several papers on user interfaces seems to indicate that long in-field studies are advised. Due to time limitations, and the very experimental design of the interface this isn't an option. The plan is to do both objective and subjective (questionnaire) measurements. Objective measurements could be e.g. time taken to complete a given task. At first I figured I needed a **paired t-test**. However, I will most likely take multiple measurements and not just one. Some reading mentioned a repeated measure ANOVA for such a scenario, but I'm not sure about this at all. Which test should I use? Which power analysis method should I use? I want to determine how many subjects I need (or can do) and what the power of my test will be.


title: Generalized Linear Model in SPSS with common values among predictors treated as subpopulations. Why?

text: I am teaching a class on logistic regression with SPSS. The textbook supplies a sample data set with a binary predictor and two numeric covariates. The sample contains 1000 rows and a number of these entries have common values for both predictors. One predictor takes only 5 values, for example, and the other takes around 20 distinct values. According to the documentation of SPSS, when this happens, SPSS treats the data as coming from subpopulations, defined through the common values. This seems to produce a different likelihood and different degrees of freedom for the AIC than what you get if you ignore subpopulations. I ran the data set in R, using glm. Degrees of freedom were 997, AIC=508.93 On SPSS, I get 99 degrees of freedom (for goodness of fit purposes) and AIC=181.341. The coefficient estimates are the same in both applications. To make matters worse, when I fit the model in SPSS with only 1 of the 2 predictors, the likelihood is LARGER than with the 2 predictor model: -87 for the 2 parameter model, and -47 for the 1 parameter model. The AIC is also dramatically smaller in the 1 parameter model, but everything else suggests that both predictors are significant and necessary. So much for the AIC criterion. I jittered the data in R, and sent it back to SPSS. I then got much the same results as in R with glm, since there were no phantom "subpopulations" for SPSS to cope with. Questions:   1. can someone supply a reference to justify treating the data as coming from subpopulations (which they actually don't in this case) when the predictors contain common value sets?   2. How am I supposed to teach model testing by comparing the deviance between two models, using SPSS and this data set, given what's going on?   3. Can I make SPSS behave like R?


title: Why is the posterior proportional to $P(X=x|C=c_i)P(C=c_i)$?

text: When we are trying to maximum a posterior, we apply Bayesian rule to convert them into posterior probabilities! $P(C=c_i|X=x) = P(X=x|C=c_i)P(C=c_i) / P(X=x)$ proportional to $P(X=x|C=c_i)P(C=c_i)$, but how can i be sure that it is really the case?


title: Analysis and writing up results

text: I am doing research for my Masters studies and am battling a bit with the statistical side of my study. I used a survey and have captured all the data into SPSS. I have no statistical background, I would thus need a bit of advice on how to interpret / analyse results. I am making use of Chi-Square test to look at relationships between categorical and other categorical variables and Correlation test for relationships between ordinal and other ordinal variables. Because of the dichotomy, I am making use of Binomial logistic regression. What I don't know is how to write results up and analyse results. For example a correlation of .000 is highly significant, but how do I word it? And what if almost all the results in the table are highly significant? For example people who agree that shopping online is safe also agree that the internet is trustworthy. Kind of expected - the result is significant but rather obvious.


title: Jointly normal distribution

text: I encountered the following sentence while reading a textbook and could not figure out why. > If $x$ and $y$ are jointly normal, then $y$ is a linear (affine) function of > $x$: $y = a+bx$. If someone could explain this to me, that would be great.


title: Confidence interval on Poisson distributed data

text: We have events which has different durations. I'm now computing confidence interval on defective events, meaning the ratio of events having encountered an issue. I consider this problem as counting defective units, so Binomial, so 1 sample proportion test. Now, I need to change my mind and compute those intervals as issues per event hour. Since some events are less than one hour, and some more, it doesn't makes sense to talk about "defective hours". So, no it is about counting defects per unit, so Poisson. Here's my question: what is the equivalent test of 1 sample proportion for Poisson distributed data? Thanks a lot for your help. P.S.: By "event", I mean a show, an exhibition, etc... Not a "statistical" event


title: Understanding the distribution of data?

text: Please help me to understand the type of the distribution of my data (width values shown below). I am providing three graphs: 1) scatter plot; 2) histogram; and 3) qqplot (variable width). Interval_number is just a serial number. Width is nothing but the time in seconds for arriving a set of 100 things. My data:               interval_number  width     1                   34      2                  234      3                  124     ....                ..      ![plot](http://i.stack.imgur.com/25ARu.png) ![histogram](http://i.stack.imgur.com/u2x03.png) ![qqnorm](http://i.stack.imgur.com/faQPa.png)


title: Inferential analysis on categorical medical data

text: I am currently conducting a retrospective chart review where I am interested in the efficacy of certain drugs when it comes to treating a certain disease (let's call the disease "A") Multiple drugs have been used to treat disease A. However, only two drugs have been so frequently used that an inferential method of analysis would generate a meaningful output. So I have two drugs, drug X and drug Y, which according to previous research are also the most efficacious drugs when treating disease A. The outcomes are coded in two dichotomous variables where 1=Success and 0=Fail. The first dichotomous variable is whether or not the drug achieved complete symptom control, and the second variable whether the symptoms were relieved significantly (defined as greater than 50% reduction of symptom frequency). So I have drug X and drug Y with outcomes coded in dichotomous variables and I would like to test the hypothesis that these drugs significantly differ in treatment outcomes (both on complete symptom control and >50% symptom reduction). I would appreciate any and all help with this problem. Thank you!


title: Suggestion for Books on Visualization and Natural Language Processing

text: ![http://nirvacana.com/thoughts/becoming-a-data- scientist/](http://i.stack.imgur.com/9f53p.png) This graphic is due to Swami Chandrasekaran. I just started learning skills for data science recently. I have basics in R, Python, Statistics, Machine Learning. Now, I want to learn the Visualization and NLP part. What books should I begin with, that covers most of what listed above? Practical books are preferred, because I want to try immediately with data available on internet (e.g. data.gov). Nevertheless, suggestions for theoretical books are welcomed too.


title: Question on Survival Analysis

text: I'm working on a survival analysis project. In part of the project, I'm asked to "identify factors that may contribute to an increased risk of death". I have difficult time understanding what I'm asked for. Could anyone help me understand this? How should I identify those factors? Is this an example of building exploratory model? If so, should I only pick variables that increase the risk? Thanks for your help.